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PREFACE 


This Annual Status Report describes the results of work performed during the 
first year of the NASA Hot Section Technology program, "3-D Inelastic Analysis 
Methods for Hot Section Components" (contract NAS3-23697). The goal of the 
program is to develop computer codes which permit more accurate and efficient 
structural analyses of gas turbine blades, vanes, and combustor liners. The 
program is being conducted under the direction of Dr. C. C. Chamis of the 
NASA-Lewis Research Center. Prime contractor activities at United Technologies 
Corporation are managed by Dr. E. S. Todd. Subcontractor efforts at the United 
Technologies Research Center, MARC Analysis Research Corporation, and the 
State University of New York at Buffalo are led by Dr. B. N. Cassenti , Dr. J. 
C. Nagtegaal , and Dr. P. K. Banerjee, respectively. 


i 



Table of Contents 


Section Page 


1.0 INTRODUCTION 1-1 

2.0 SUMMARY 2-1 

2.1 Constitutive Models 2-1 

2.2 Mechanics of Materials Model 2-2 

2.3 Special Finite Element Model 2-2 

2.4 Advanced Formulation (Boundary Element) Model 2-2 

3.0 TECHNICAL PROGRESS 3.1-1 

3.1 Constitutive Models 3.1-1 

3.1.1 Simplified Secant Elastic Model 3.1-1 

3.1.2 Current State-of-the-Art Model 3.1-5 

3.1.3 Advanced Viscoplastic Model 3.1-12 

3.1.4 List of Symbols 3.1-15 

3.1.5 References 3.1-17 

3.2 Mechanics of Materials Model 3.2-1 

3.2.1 Computer Program: Formulation/Description 3.2-1 

3.2.2 Program Validation/Verification 3.2-5 

3.2.3 List of Input Parameters 3.2-10 

3.2.4 List of Symbols 3.2-13 

3.3 Special Finite Element Model 3.3-1 

3.3.1 Literature Survey 3,3-1 

3.3.2 Formulation Development 3,3-3 

3. 3. 2.1 Introduction 3.3-3 

3. 3. 2. 2 Mixed Forms and Iterative Solutions 3.3-4 

3. 3. 2. 3 Element Technology 3.3-8 

3. 3. 2. 4 Inelastic Constitutive Models 3.3-10 

3. 3. 2. 5 Time Integration 3.3-11 

3. 3. 2. 6 Eigenvalue Extraction 3.3-12 

3.3.3 Program Development 3.3-12 

3. 3. 3.1 Introduction 3.3-12 

3. 3. 3. 2 Overview and Control Structure 3.3-13 

3. 3.3. 3 Element Library 3.3-17 

3. 3. 3. 4 Nonlinear Analysis Capabilities 3.3-17 

3.3.4 Program Validation/Verification 3.3-18 

3.3.5 Input Data Structure for the MHOST Program 3.3-33 

3.3.6 List of Symbols 3.3-37 

3.3.7 References 3.3-39 


ii 



Table of Contents (continued) 

Secti on Page 


3.4 Boundary Element Method 3.4-1 

3.4.1 Overview 3.4-1 

3.4.2 Literature Survey 3.4-1 

3. 4. 2.1 Introduction 3.4-1 

3. 4. 2. 2 Linear Stress Analysis 3.4-2 

3.4. 2. 3 Dynamic Stress Analysis 3.4-3 

3. 4. 2. 4 Nonlinear Stress Analysis 3.4-4 

3.4.3 Formulation Development 3.4-6 

3. 4. 3.1 Summary 3.4-6 

3. 4. 3. 2 Linear and Nonlinear Stress Analysis 3.4-6 

3. 4. 3. 3 Dynamic Stress Analysis 3.4-13 

3.4.4 Computer Program Development 3.4-19 

3.4. 4.1 Introduction 3.4-19 

3. 4. 4.2 Global Program Structure 3.4-19 

3. 4. 4. 3 Program Input 3,4-22 

3. 4. 4. 4 Surface Integral Calculation 3.4-25 

3. 4. 4. 5 Volume Integral Calculations 3.4-28 

3. 4. 4. 6 System Matrix Assembly 3.4-29 

3. 4. 4. 7 System Equation Solution 3,4-31 

3. 4. 4. 8 Inelastic Solution Process 3.4-31 

3.4.4. 9 Output 3.4-32 

3.4.5 Program Validation/Verification 3,4-33 

3.4.5. 1 Validation of Elastic Capabilities 3.4-33 

3.4. 5. 2 Validation of Inelastic Algorithms 3.4-39 

3. 4. 5. 3 Validation of Dynamic Analysis 3,4-47 

3. 4. 5. 4 Notched Specimen Verification 3.4-50 

3. 4. 5. 5 Hot Section Component Analysis 3.4-56 

3.4.6 Boundary Element Stress Technology (BEST) 3,4-59 

Program U911- Input Description 

3. 4. 6.1 Case Control Input 3.4-60 

3. 4. 6. 2 Material Property Input 3.4-60 

3.4. 6.3 Generic Modeling Region Input 3.4-62 

3. 4. 6. 4 Interface Input 3.4-63 

3.4. 6. 5 Boundary Condition Set Input 3.4-64 

3. 4. 6. 6 Body Force Input 3.4-65 

3.4.7 Sample Output from BEST 3.4-66 

3. 4. 7.1 Input Echo 3,4-66 

3. 4. 7. 2 Case Control Summary 3.4-67 

3. 4. 7. 3 Generic Modeling Region (GMR) Definition 3.4-67 

3.4. 7.4 Boundary Condition Definition 3.4-68 

3. 4. 7. 5 Boundary Solution (Element Basis) 3.4-68 

3.4. 7.6 Boundary Solution (Nodal Basis) 3.4-69 

3.4.7. 7 Cell Node Displacements 3.4-70 

3. 4. 7. 8 Cell Node Stresses 3.4-70 

3.4. 7. 9 Cell Node Strains 3.4-71 

3.4.8 List of Symbols 3.4-72 

3.4.9 References 3.4-74 


DISTRIBUTION LIST 



SECTION 1.0 


INTRODUCTION 


Aircraft povverplant fuel consumption and expenditures for repair/replacement 
of worn or damaged parts make up a significant portion of commercial avia- 
tion's direct operating costs. For modern gas turbines, both factors depend 
heavily on the degree to which elevated flowpath temperatures are sustained in 
the hot section modules of the engine. Higher temperatures reduce fuel con- 
sumption by raising the basic efficiency of the gas generator thermodynamic 
cycle. At the same time, these elevated temperatures work to degrade the dura- 
bility of structural components (combustor liners, turbine blades and vanes, 
airseals, etc.) that must function adjacent to or within the hot gaspath it- 
self, leading in turn to larger maintenance/material costs. Pursuit of the 
best compromise between performance and durability presents a challenge that 
will continue to tax the ingenuity of advanced gas turbine design analysts for 
years to come. 

Hot section durability problems appear in a variety of forms, ranging from 
oxidation/corrosion, erosion, and distortion (creep deformations) to occur- 
rence of fatigue cracking. Even modest changes in shape, from erosion or dis- 
tortion of airfoils for example, can lead to measurable performance deteriora- 
tion that must be accurately predicted during propulsion system design to in- 
sure that long-term efficiency guarantees can be met. Larger distortions in- 
troduce serious problems such as hot spots and profile shifts resulting from 
diversion of cooling air, high vibratory stresses associated with loose tur- 
bine blade shrouds, difficult disassembly/reassembly of mating parts at over- 
haul, etc. These problems must be considered and efforts made to eliminate 
their effects during the engine design/development process. Initiation and 
propagation of fatigue cracks represents a direct threat to component struc- 
tural integrity and must be thoroughly understood and accurately predicted to 
insure continued safe and efficient engine operation. 

Accurate prediction of component fatigue lives is strongly dependent on the 
success with which inelastic stress/strain states in the vicinity of holes, 
fillets, welds, and other discontinuities can be calculated. Stress/strain 
computations for hot section components are made particularly difficult by two 
factors - the high degree of geometrical irregularity which accompanies so- 
phisticated cooling schemes, and complex nonlinear material behavior associ- 
ated with high temperature creep/plasticity effects. Since cooling air ex- 
traction reduces engine cycle efficiency, concerted efforts are made to mini- 
mize its use with the result that elaborate internal passages and surface 
ports are employed to selectively bathe local regions (airfoil leading edges, 
louver liner lips, etc.) for which the high temperature environment is most 
severe. These cooling features frequently interrupt load paths and introduce 
complex temperature gradients to the extent that the basic assumptions of one- 
and two-dimensional stress analysis procedures are seriously compromised and 
the use of three-dimensional techniques becomes mandatory. Even in the pres- 
ence of cooling, component temperature and stress levels remain high relative 
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to the material's melting point and yield strength values. The combinations of 
centrifugal, aerodynamic, thermal, and other mechanical loadings that typical- 
ly occur in flight operation then serve to drive the underlying material re- 
sponse beyond accepted limits for linear elastic behavior and into the regime 
characterized by inelastic, time-dependent structural deformations. Thus, an 
ability to account for both complexities, three-dimensional and inelastic 
effects, becomes essential to the design of durable hot section components. 

General purpose finite element computer codes containing a variety of three- 
dimensional (brick) elements and inelastic material models have been available 
for more than a decade. Incorporation of such codes into the hot section de- 
sign process has been severely limited by high costs associated with the ex- 
tensive labor/computer/time resources required to obtain reasonably detailed 
results. Geometric modeling systems and automated input/output data processing 
packages have received first attention from software developers in recent 
years and will soon mature to the point that previous over-riding manpower 
concerns will be alleviated. Prohibitive amounts of Central Processing Unit 
(CPU) time are still required for execution of even modest-size three-dimen- 
sional inelastic stress analyses, however, and is chief among the obstacles 
remaining to be remedied. With today's computers and solution algorithms, 
models described by a few hundred displacement degrees of freedom commonly 
consume one to three hours of mainframe CPU time during simulation of a single 
thermomechanical loading cycle. A sequence of many such cycles may, of course, 
be needed to reach the stabilized conditions of interest. Since accurate 
idealizations of components with only a few geometrical discontinuities can 
easily contain several thousand degrees of freedom, inelastic analysis of hot 
section hardv/are with existing codes falls outside the realm of practicality. 

The Inelastic Methods program addresses the need to develop more efficient and 
accurate three-dimensional inelastic structural analysis procedures for gas 
turbine hot section components. A series of new, increasingly rigorous, stand- 
alone computer codes is being created for the comprehensive numerical analysis 
of combustor liners, turbine blades and vanes. Theoretical foundations for the 
codes feature mechanics of materials models, special finite element models, 
and boundary element models. Heavy attention will be given to evolution of 
novel modeling methods that permit non-burdensome yet accurate representations 
of geometrical discontinuities such as cooling holes and coating cracks. A 
selection of constitutive relations has been provided for economical or so- 
phisticated description of inelastic material behavior as desired. Finally, 
advantages which accrue from application of the improved codes to actual com- 
ponents will be demonstrated by execution of benchmark analyses for which 
experimental data exist. 
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SECTION 2.0 


SUMMARY 


The 3-D Inelastic Analysis Methods program is divided into two 24-month seg- 
ments: a base program, and an option program to be exercised at the discretion 
of the Government. During the base program, a series of new computer codes em- 
bodying a progression of mathematical models (mechanics of materials, special 
finite element, boundary element) is being developed for the streamlined anal- 
ysis of combustor liners, turbine blades and turbine vanes. These models will 
address the effects of high temperatures and thermal /mechanical loadings on 
the local (stress/strain) and global (dynamics, buckling) structural behavior 
of the three selected components. 

The first year (Task I) of the base program dealt with "linear" theory in the 
sense that stresses/strains and temperatures in generic modeling regions are 
linear functions of the spatial coordinates, and solution increments for load, 
temperature and/or time are extrapolated linearly from previous information. 
Three linear formulation computer codes, hereafter referred to as MOMM (Me- 
chanics £f Materials Model), MHOST (MARC -HOST ), and BEST (^oundary Hement 
Stress Technology), have been created and are described in more detail in 
"Sections 3.2, 3.3, and 3.4, respectively. 

The second half of the base program (Task II), as well as the option program 
(Tasks IV and V), will extend the models to include higher-order representa- 
tions of deformations and loads in space and time to deal more effectively 
with collections of discontinuities such as cooling holes and coating cracks. 
Work on Task II (polynomial theory) has commenced and will be the subject of 
primary interest in the next Annual Status Report. 

2.1 CONSTITUTIVE MODELS 

Three increasingly rigorous constitutive relationships are employed by MOMM, 
MHOST, and BEST to account for nonlinear material behavior (creep/plasticity 
effects) in the elevated temperature regime. The simplified model assumes a 
bilinear approximation of stress-strain response and generally glosses over 
the complications associated with strain rate effects, etc. (Section 3.1.1). 
The state-of-the-art model partitions time-independent (plasticity) and time- 
dependent (creep) phenomena in the conventional way, invoking the Mises yield 
criterion and standard (isotropic, kinematic, combined) hardening rules for 
the former and a power law for the latter (Section 3.1.2). Walker's viscoplas- 
tic theory, which accounts for the interaction between creep and plasticity 
that occurs under cyclic loading conditions, has been adopted as the advanced 
constitutive model (Section 3.1.3). 
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2.2 MECHANICS OF MATERIALS MODEL 

In essence, the Mechanics £f Materials Model (MOMM) is a stiffness method fi- 
nite element code that utilizes one-, two- and three-dimensional arrays of 
beam elements to characterize hot section component behavior. Limitations of 
such beam model representations are recognized, of course, but are fully ac- 
ceptable in view of the benefits of having a fast, easy to use, computation- 
ally efficient tool available for application during the early phases of com- 
ponent design. The full complement of structural analysis types (static, buck- 
ling,, vibration, dynamics) is provided by MOMM, in conjunction with the three 
constitutive models mentioned above. Capabilities of the code have been tested 
for a variety of relatively simple problem discretizations (examples are pro- 
vided in Section 3.2.2). Work to establish modeling guidelines for simulation 
of two- and three-dimensional behavior is in progress. 

2.3 SPECIAL FINITE ELEMENT MODEL 

The MHOST (MARC -HOST ) code employs both shell and solid (brick) elements in a 
mixed method framework to provide comprehensive capabilities for investigating 
local (stress/strain) and global (vibration, buckling) behavior of hot section 
components. Development of the code has taken full advantage of the wealth of 
technical expertise accumulated at the MARC Corporation over the last decade 
in support of their own commercially available software packages to create 
new/improved algorithms (Section 3.3.2) that promise to significantly reduce 
CPU (central processing unit) time requirements for three-dimensional analy- 
ses. First generation (Task I) MHOST code is operational and has been tested 
with a variety of academic as well as engine-related confiaurations (Section 
3.3.4). 

2.4 ADVANCED FORMULATION (BOUNDARY ELEMENT) MODEL 

Successful assembly of the all -new BEST (£oundary Hement £tress Technology) 
code constitutes perhaps the most important accomplishment of the Task I ef- 
fort. The difficult challenge of extending the basic theory and algorithms to 
encompass inelastic and dynamic effects in three-space was effectively met by 
combining the special skills and efforts of the research and programming teams 
at SUNY-B and P&W. As with MOMM and MHOST, the initial version of BEST is exe- 
cutable and has been exercised with a number of small and large test cases 
(Section 3.4.5). While MHOST and BEST are currently viewed as mutually com- 
plementary, they are also competitors; and overall performance on large in- 
elastic models will be watched with high interest as the codes mature. 
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SECTION 3,0 


TECHNICAL PROGRESS 


3,1 CONSTITUTIVE MODELS 

Three material models are available for use with the mechanics of materials, 
special finite element, and boundary element models: 1) a simplified material 
model, 2) a state-of-the-art material model, and 3) an advanced material mod- 
el, The simplified model uses secant moduli and assumes a bilinear stress- 
strain response which is currently neither strain-rate nor temperature depen- 
dent, Later versions of the simplified material model will include provisions 
for both temperature and strain-rate dependence. The state-of-the-art material 
model is a standard elastic-plastic-creep model (Reference 1), The advanced 
model is a modified form of Walker's viscoplastic material model (References 2 
and 3), The following sections provide a detailed discussion of each of these 
models. 


3,1,1 Simplified Secant Elastic Model 


In the simplified elastic model, stress-strain curves for various strain rates 
are the basic input material properties. Tension response is assumed to be the 
same as compression response. The initial response is represented by an elas- 
tic material with modulus, Eq, and Poisson's ratio, vq. At the conclusion 
of the calculation for the response, an equivalent strain is predicted. At 
this strain, two equivalent stresses can be considered: 1) the calculated 
stress, and 2) the stress from the input stress-strain curves at the predicted 
strain. If the two stresses are sufficiently close in value, then the calcula- 
tions can be terminated. If the two stresses are not sufficiently close, then 
the new modulus is taken to be the stress from the stress-strain curves divid- 
ed by the strain, and the calculations are repeated. 

This concept must now be expanded to multidimensional stress states. For this 
purpose, consider an elastic material, then: 
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The stress and strain can be partitioned into deviatoric and volumetric parts. 


^ij = ®ij ^kk S'j 

(3.1-2) 

°ij = ^ij “"kk 

(3.1-3) 


The volumetric components, from equation (3.1-1) are related by 
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where K is the bulk modulus. 

The deviatoric parts can be shown to be related by 
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where J2 ‘'s the second invariant of the deviatoric strain tensor. 
Equation (3.1-7) now becomes 
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Since only the ratio a/F will be used to represent the material response, an 
additional assumption is needed to obtain the second elastic constant. For 
this purpose, assume the bulk modulus is constant, and given by equation 
(3.1-4) 
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TTr-"2vT 




2 (1 ^ v^) 6^ 
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(3.1-10) 


where Gq, Eq, vq are the moduli and Poisson's ratio at the origin (i.e., 
F=7=0). The current shear modulus is known from the slope F/F. Then from equa- 
tion (3.1-9) 
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Solving equation (3.1-11) for 
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Figure 3.1-1 presents the variation in Poisson's ratio with modulus. The 
Young's modulus can be determined from 


E = 2 (1 + v) G = (1 + v) a/e 


As an example, consider a uniaxial stress state 
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o = a and e = (1 + v) E = Ej^j^ - £22 (3.1-15) 

The F, F curve is now the input stress-strain curve. 

To illustrate the convergence of the iterative procedure, consider three par- 
allel bars supporting an equivalent total load. The bars are assumed to be 
elastic-plastic. Each has a Young's modulus of 10 x 10° psi and a hardening 
slope of 0.5 X 10® psi. The yield stresses are different. The central bar 
will be assumed to have a yield stress of 20 ksi while the two outer bars have 
a yield stress of 10 ksi. The area of each bar is 1/3 in^, making a total 
area of 1.0 in^. Figure 3.1-2 illustrates that convergence has occurred in 
six iterations for a total load of 30,000 lb and that each of the bars has 
yielded. 

The material constants for the simplified model are input to the computer code 
through data input cards. 
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3.1.2 Current State-of-the-Art Model 

The current state-of-the-art model has been taken to be the classical elastic- 
plastic-creep model that is available in the MARC code, and described in Ref- 
erence 1. The creep model is essentially a steady state power law (stress) 
model. The plasticity model includes isotropic, kinematic, and a combined 
hardening law. Both the creep and plasticity models assume no permanent volu- 
metric deformations. For the mechanics of materials computer code, the materi- 
al properties for the state-of-the-art constitutive model are included in data 
statements in subroutine SOACON. 

Plastic Iteration Procedure 


Consider the case of a small strain elastic-plastic response of a typical 
structure. Sufficiently large applied loads will result in permanent or plas- 
tic deformation. A procedure for calculating the response of the structure un- 
dergoing plastic deformation is required. 

To evaluate the response of the structure, the loading history is divided into 
a number of incrementally applied loading steps. Each of these load increments 
can then be applied sequentially to the structure. An iterative scheme is then 
required to calculate the response of the structure to each individual load 
increment. 

At the beginning of a new load increment it may be assumed that the strain 
will change in a manner analogous to the previous increment. As an initial 
estimate all of the strain change is then assumed to be elastic. The change in 
the stresses can then be calculated using Hooke's Law or 


^‘"ij “ '-ijkl^^kl 


where: Aa.. 

* V 


is the incremental stress vector. 


A€|^^ is the incremental total strain vector, and 


(3.1-16) 


L 


e 

ijkl 


is the matrix of elastic constants. 


If the resulting total stress is within the yield surface, the matrix of mate- 
rial constants, L-jji^] is simply given by 



(3.1-17) 


If the resulting total stress is outside the yield surface, weighted material 
constants and stiffness matrices will have to be calculated. It should be 
noted at this point that if a load increment is exceedingly large and if there 
is a sudden change in the type of loading, care must be taken in order to 
iterate to the correct solution. 
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If the resulting total stress is outside the yield surface, the fraction of 
the stress increment that remains elastic must be determined. This corresponds 
el 

to Ae. .l in Figure 3.1-3. I,f the yield surface in stress space is considered 

* %I 

to be given by 


«»ij> 


0 . 


then the appropriate m in 

f(aj:^ + maa..) = 0 (3.1-18) 

I J I J 

may be determined where is the stress tensor from the previous increment. 

The mean material matrix is calculated from 


^ijkl = "^ijkl ‘■ijkl (3.1-19) 

where is the tensor relating and . 

Once the tensor has been determined, standard solutions can be applied 

to find the incremental changes in the displacements, strains and loads. For 
example, if the strains are given by 

|Ae} = [b] I&u} (3.1-20) 

where |au} is the vector of incremental nodal displacements, and Cb] is the 
matrix relating the vector of element strains jaef to the nodal displacements, 
the stiffness matrix can be found from 

CK] = f [b]"^ CD] Cs] dV (3.1-21) 

V 

where [D] is the matrix representation of the tensor Lfjid. The strain-dis- 
placement matrix [b] depends on the formulation of the problem. 

The incremental nodal displacements and strains can be evaluated by solving 
for Au in 


[K] ]au| = |aP[ + {aG} (3.1-22) 

and then applying equation (3,1-20). 

In the mechanics of materials computer code the stiffness matrix K is held 
constant, and changes in the stiffness matrix are included in aG. 
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Figure 3,1-3 Elasti'c-Plastic Strain Decompositions for Bilinear Stress-Strain 
Law 


The term aP in equation (3,1-22) is the applied incremental load. The term aG 
is defined as the pseudo-load correction to the stiffness matrix due to in- 
elastic strains which is added to equation (3.1-22). The aG vector calculated 
from creep strain, for example, is shown in equation (3.1-34). 

One iteration cycle is completed each time the stiffness matrix is formed and 
the resulting equations solved. At the end of each cycle the resulting solu- 
tion must be tested for convergence. This is accomplished, by considering the 
change in energy. 


r 



gN _ ^N-1 

1/2 (E^ + 


(3.1-23) 


where E*^"^ is the change in energy summed over all elements on the previous 
cycle and E^ is the energy including the present cycle. 


An accurate solution will usually result if r is maintained less than 0.1 for 
elastic-plastic problems. 
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If the solution has satisfied the convergence, the stresses and strains can be 
updated and a new load increment added. If the solution has not converged, 
then a new guess for the strains, based on the latest cycle, must be input and 
the calculation procedure repeated. When the solution has not converged after 
a given number of cycles, the program should exit from the load incrementing 
loop. 


Figure 3.1-4 is a flow chart illustrating the small strain elastic-plastic it- 
eration procedure. 


For isotropic materials the moduli in equation (3.1-19) are given by 


^ijkl = IS'k^jl'" T -'Zv S’j ^kl 


and 


(3.1-24) 


L! 


e-p 


ijkl “ 1+v l^ik^jl ^ l-2v ®ij°kl 

where E is Young's modulus 

V is Poisson's ratio 

6.. is the Kronecker delta 

* J 

. « .P 

“ij - ® 'fj 


«-• .-Si. 


3/2 (S.j -S2^.j) (S|^^ 


= 2 /I + v\ "/I + v\ 

1 + I 


(3.1-25) 


(3.1-26) 

(3.1-27) 


P 

are the plastic strain rates 


f = yJz/3 


P .P 

. e « • 

j ij 


G is the kinematic hardening slope 


(3.1-28) 


H is the isotropic hardening slope 
ay is the initial yield stress, and 
^ij = ^ij ■ 1/3 ak|<6-jj is the deviatoric stress. 


The strain rate has been decomposed into elastic (including thermal), plastic 
and creep components, or 


e 

e 


.P 

e 


+ 


.C 

e 


(3.1-29) 
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The plastic yield surface was assumed to satisfy an equivalent Mises yield 
surface given by 

1/2 (S.j. -S2.J.) (S.j -fi.j) = 1/3 djj (3.1-30) 

The method presented in Reference 4 is used to calculate the elastic-plastic 
modul i . 



Figure 3.1-4 Elastic-Plastic Iteration Procedure 
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Time Effects Iteration Procedure 


The creep strain rate will depend in general on the stress, the accumulated 
creep strain, the temperature and time. To illustrate the incrementing proce- 
dure, assume that the creep strain rate is normal to the Mises yield surface 
in stress space, then the creep strain rate is given by 


cr 


V 3/2 S 




kl -^kl 
K 


3/2 S. . 


n/3/2 S, 


S 

mn mn 


(3.1-31) 


For a specific time increment the incremental creep strains were approximated 

by 


AeT; = At. (3,1-32) 

I J I J 

The incremental displacements are 

[K] |au| = |aP[ + |aG[ (3.1-33) 

where 

|aG[ = f Cb]^ [E] Iae^} dV (3.1-34) 

is the pseudo-creep load, |AeC[ is the vector of element creep strain, and 
[E] is the elasticity matrix. The strain increment can be calculated from 
equation (3.1-20) and the strains, creep strains, stresses and displacements 
can be updated. 

A convergence test on the stresses should be performed. If the algorithm has 
not converged, a shorter time step should be used and the calculations re- 
peated. If the criterion has been satisfied, then the time step can be in- 
creased. Figure 3.1-5 is a flow chart illustrating the small strain creep it- 
eration procedure. 
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Figure 3.1-5 Creep Iteration Procedure 
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3.1.3 Advanced Viscoplastic Model 


The viscoplastic model described in References 2 and 3 has been selected as 
the advanced constitutive model. Reference 2 describes the basic theory; while 
Reference 3 describes modifications to the form of the basic theory, and modi- 
fications to the material parameters for Hastelloy X. The modifications pro- 
vide more accuracy at relatively low temperatures. 

For uniaxial loading the viscoplastic material model (Figure 3.1-6) reduces to 


sgn(. -n) * -a) H - k) <»> 

0,0 - k0 

(3.1-35) 

• « • 

£2 = ngC - n^ |c |£2 

(3.1-36) 

C = E - ^ 

(3.1-37) 


where C is the inelastic strain, 

£2 is the back stress, 
a is the stress, 
e is the strain, and 

k, CToo» n, K, n 2 , n 3 and E are material constants. 

The absolute value and unit ramp functions are represented by 

|X| ={-;;*<§ (3.1-38) 

and 

<x>={° (3.1-39) 

The inelastic strain in equation (3.1-35) consists of two components: 1) a 
time dependent power law creep component, containing the material constants n 
and k, and 2) a time independent plastic component, containing the material 
constants 0^0 and k. The parameter becomes equivalent to the yield stress 
as: 1) k, in equation (3.1-35), approaches unity, and 2) the back stress, fi, 
approaches zero. The back stress is a key variable in many viscoplastic ma- 
terial models. Its evolution is given by equation (3.1-36). Equation (3.1-37) 
represents the inelastic strain as the difference between the total strain and 
the elastic strain. 
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Figure 3.1-6 Modified Walker's Theory 


Subroutine HYPELA in the mechanics of materials computer code integrates 
Walker's viscoplastic equations and calls subroutine HYPCON to evaluate the 
material parameters. HYPCON contains the latest estimates for the parameters 
in the modified Walker's theory. Each load increment in the analysis is divid- 
ed into NSPLIT subincrements. The integration of the constitutive equations is 
performed by using forward differences with a step size determined by dividing 
the load increment by NSPLIT. Subroutine HYPELA performs the integration in 
two ways: 1) a fixed step size, or 2) a variable step size. In the fixed step 
size, forward difference NSPLIT is the same for all load increments and sub- 
increments. 
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In the variable step size, forward difference NSPLIT is determined by the mag- 
nitude of the change in a strain measure for every subincrement. The change in 
the strain measure is defined as 


3aJ» 

E = aR + (3.1-40) 

where 


aR = V 2/3 aC.. aC.. (3.1-41) 

I O * J 

aJ« = 3/2 AS.. AS., and (3.1-42) 

t ' J * J 

the quantity ae is calculated and is stored as variable ERRORO. There are 
three possible ways to determine NSPLIT. The method depends on the size of 
ERRORO. If 


ERR0R2 < ERRORO < ERRORl, (3.1-43) 

then NSPLIT remains the same for the next subincrement (ERRORl and ERR0R2 are 
user-specified in HYPELA). If 


ERRORO < ERR0R2, (3.1-44) 

then NSPLIT is divided in two for the next subincrement and rounded (up) to 
the nearest integer. If 

ERRORO > ERRORl, (3.1-45) 

then NSPLIT is doubled and the step is recomputed. The value of NSPLIT at the 
end of the increment is stored in the state variable TEMP (16). The initial 
value of NSPLIT is user-specified in HYPELA. The maximum value of NSPLIT is 
specified by MXSPLT. If NSPLIT exceeds MSXPLT, the message: 

"UNABLE TO REDUCE ERROR IN LESS THAN MXSPLT SUBINCREMENTS" 

is written where the value of MXSPLT is inserted in the WRITE statement. After 
this, the integration is performed using a constant step size. 
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3.1.4 List of Symbols 


List of Symbols 
Referenced Within Section 3.1 


Symbol 

«ij 

V 

E 

K 

J2 

J2 


Description 


Page 

Strain 


3.1-1 

Stress 


3.1-1 

Kronecker delta 


3.1-1 

Poisson's ratio 


3.1-1 

Young's modulus 


3.1-1 

Deviatoric strain 


3.1-2 

Deviatoric stress 


3.1-2 

Bulk modulus 


3.1-2 

Second invariant of the 
stress tensor 

deviatoric 

3.1-2 

Second invariant of the 
strain tensor 

deviatoric 

3.1-2 


Equivalent stress 


3.1-2 


e 

Equivalent strain 

G 

Shear modulus 

1-ijkJj 

Matrix of material constants 

|au| 

Incremental nodal displacements 

]ac[ 

Incremental stress 

[b] 

Strain- displacement matrix 

1 — 1 
I 1 

Stiffness matrix 
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3.1- 2 

3.1- 5 

3.1- 6 

3.1- 5 

3.1- 6 
3.1-6 
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List of Symbols 
Referenced Within Section 3.1 


Symbol 

Description 

Page 

|ap} 

Incremental applied load vector 

3.1-6 

|aG} 

Incremental pseudo-load vector 

3.1-6 

en 

Energy in Nth cycle 

3.1-7 

r 

Convergence parameter 

3.1-7 

CE] 

Elasticity matrix 

3.1-10 

G 

Kinematic hardening slope 

3.1-8 

H 

Isotropic hardening slope 

3,1-8 

S2 

Back stress 

3.1-12 

c 

Inelastic strain 

3.1-12 

0 

X j M > n, \ 
m, ni, 02 , / 
03. 14, ns, > 

Material constants 

3.1-12 


K2 ) k » O'oo 
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3.2 MECHANICS OF MATERIALS MODEL 


3.2.1 Computer Program: Formulation/Description 

The three-dimensional nonlinear mechanics of materials finite element computer 
program utilizes an intersecting network of beams to model a structural com- 
ponent. The program calculates the total strain as a linear function of posi- 
tion in cross section and along the length of the beam. Three material models 
are included in the code: the simplified material model. Walker's viscoplastic 
material model, and the state-of-the-art material model. Static and transient 
analyses can be performed with applied loads, thermal loads, and enforced dis- 
placements. The lowest frequency and mode shape using either initial or tan- 
gent stiffness is calculated; and buckling analysis is included in the static 
problem using initial or tangent stiffness. The program flow is summarized in 
Figure 3.2-1. 

Input parameters to the computer code consist of information defining the mod- 
el itself and information describing the method of solution desired. The model 
is defined by beams which are connected at grid points. The element coordinate 
system of a given beam is defined by an orientation grid point. The geometry 
of a beam is rectangular in cross section, with the dimensions of the cross 
section along the element coordinate axes specified. The material properties 
are specified for each beam, including Young's modulus, Poisson's ratio, mass 
density, coefficient of expansion, and yield stress. The initial temperature 
of the beam network is input, and the time at initial conditions is set to 
zero. A hardening slope for use with the simplified material model is entered, 
with a zero slope indicating perfectly-plastic behavior. Boundary conditions 
are specified by indicating at each node, excluding orientation grid points, a 
constrained or nonconstrained condition for the six degrees of freedom. 

Input associated with the selection of the method of solution include the pa- 
rameters that indicate: 

1. the choice of constitutive model to be used, 

2. the choice of a static or transient analysis, 

3. the choice of initial or tangent stiffness in solving for the lowest 
frequency and mode shape, and 

4. the choice of including buckling analysis with either initial or tan- 
gent stiffness. 

The number of integration points in each beam is user-specified; stresses and 
strains are calculated at each integration point, and the user specifies be- 
tween two and ten points along each element coordinate axis direction in each 
beam. The convergence value, defining the allowable energy change between two 
consecutive iterations in the static analysis or allowable range in internal 
energy for the adaptive time step calculation in the transient analysis, is 
entered by the user. The number and Jtype of loading increments are also speci- 
fied. 


3.2-1 




Figure 3.2-1 
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The stiffness and mass matrices for each beam in the element coordinate system 
are computed and transformed to the global coordinate system. The stiffness 
and mass matrices are then assembled to form the global mass and stiffness ma- 
trices. The boundary conditions are applied to the stiffness matrix, and the 
matrix is then inverted. Any change in the stiffness due to nonlinear effects 
will be accounted for in the pseudo-load vector; therefore, the stiffness ma- 
trix is only inverted once. 

Depending on user- input, the program now is directed to the appropriate branch 
of the program: static or transient analysis. For static analysis, the loading 
increment is read from the data input, including forces and moments or en- 
forced displacements, specified at each degree of freedom of the structure. 
The temperature increment is also entered. An initial incremental displacement 
vector is set to zero and strain, stress and pseudo-load vectors are calculat- 
ed from the incremental displacement vector using the mechanics of materials 
model selected by the user. The pseudo-load vector accounts for the effects of 
nonlinearity and allows the use of the original stiffness matrix throughout 
the calculations. The equations governing the system are as follows: 

[K] |au} = |aP[ + |aG[ (3,2-1) 

where [K] = elastic stiffness matrix, 

|au| = incremental displacement vector, 

laPf = incremental applied load vector, and 

|aG} = pseudo-load vector, due to inelastic strains. 


|aG} = /[ef [E] |ae} dV (3.2-2) 

where [B] = strain-displacement matrix, 

[E] = elasticity matrix, and 

|ae} = inelastic strains. 

Equation (3.2-1) is solved for the incremental displacement vector, au, which 
is substituted for the initial incremental displacement vector and used in the 
second iteration, continuing until the change in energy in two consecutive it- 
erations is less than the convergence value input by the user. When conver- 
gence occurs, the incremental loading, displacements, strains and stresses for 
that loading increment are printed; the total load, displacement, strain and 
stress vectors, as well as temperature, are then updated. Each loading incre- 
ment is read in and executed similarly, and the values of stress, strain and 
displacement for the total loading are calculated and printed upon conclusion 
of the last increment. 

The transient analysis is based on a simple Euler integration and includes a 
self-adaptive time step scheme. Damping is not included directly in the tran- 
sient analysis but is present in the viscoplastic material models. The loading 
for each increment is the total load at that given time, which is entered into 
the program by a user-supplied subroutine. The temperature increment and time 
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step are also entered. As in the static branch, the initial displacement vec- 
tor is set to zero and the strains, stresses and pseudo-load vector are cal- 
culated using the designated mechanics of material model. An Euler integration 
is then used to calculate current displacements at the end of the present time 
step. The governing equations are as follows: 



|a} = ]aF} - [K] IaU^} 

(3.2-3) 


]aV} = [M]’^ {a} * DT 

(3.2-4) 


IaUjI = (|V| 1/2 |AVf) * DT 

(3.2-5) 

where |a[ 

= acceleration vector. 


|aF[ 

= applied and pseudo-loads. 


[l(] 

= elastic stiffness matrix. 



|aUq} = displacement vector at beginning of time step, 
|aV[ = change-in-velocity vector. 


[M] = mass matrix, 

DT = time step, 

|aUj[ = displacement vector at end of time step, and 
|v| = velocity vector. 

A measure of the work done and the change in internal energy of the system 
during the time step is computed, and the time step is adapted accordingly. If 
the time step is accepted, the current displacements, strains and stresses are 
printed, and the current displacements are inserted for the initial displace- 
ments in the following time step. If the time step is unacceptable according 
to the adaptive scheme, the time step is changed, the load is recalculated, 
and the displacements are reset to the initial value at the beginning of that 
time step. The analysis continues until the user-designated number of incre- 
ments is completed. 

Following the static or transient analysis, the user has a choice between cal- 
culating the lowest frequency and mode shape or all frequencies and mode 
shapes. The method of solution for the calculation of the lowest frequency and 
mode shape is the inverse power method, which is represented by the following 
expression: 

(CK]-^ CM] - xCl]) {x.^J = |x^| (3.2-6) 

where [K] = stiffness matrix, 

[M] = mass matrix, 

[I] = identity matrix, 

X = eigenvalue, and 
= eigenvector. 

The method of solution in the calculation of all frequencies and mode shapes 
for a given problem is the Jacobi method, which is based on simple similarity 
transformations. 
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The procedure for determining the coefficients of the inverse stiffness matrix 
is one that can represent the original stiffness of the structure or the cur- 
rent stiffness including nonlinear effects. A small load is placed at one of 
the nonconstrained degrees of freedom of the structure, and the displacements 
are computed using the specified constitutive model. The coefficients of the 
appropriate row of the inverse stiffness matrix are calculated by dividing the 
calculated displacements by the applied force. This procedure is continued for 
each nonconstrained degree of freedom until an inverse stiffness matrix, with 
dimensions equal to the number of nonconstrained degrees of freedom of the 
structure, is formed. If the frequency is to be calculated using the initial 
stiffness of the structure, all variables used in the static or transient 
analysis are set to the original values. If the tangent stiffness is request- 
ed, all variables retain the current values for use in the frequency calcula- 
tion. Only the initial stiffness option is available for use in a transient 
analysis since current stiffness cannot be readily calculated. 

Buckling analysis can be executed in a static problem. The buckling analysis 
is based on a two step process similar to that in the NASTRAN finite element 
code. In the first step, the beam loads are determined. In the second step, a 
first order large displacement correction, proportional to the loads, is in- 
cluded in the stiffness matrix. Buckling occurs when the determinant of the 
new matrix vanishes. In the determination of the stiffness matrix used in the 
buckling calculation, the stiffness coefficients are calculated in the same 
fashion as was described in the frequency calculation, with the user choosing 
the initial or tangent stiffness. The beam loads are calculated using the ini- 
tial stiffness matrix and then adding the pseudo-load vector. The actual buck- 
ling calculations are accomplished using the inverse power method to find the 
critical buckling factor and the buckled shape. 

3.2.2 Program Validation/Verification 

Some of the test cases (i.e., TESTl - TEST5) which have been executed to vali- 
date the Mechanics o^f Materials Model (MOMM) computer code are summarized be- 
low. Each of these cases test varTous segments of the theory and computer code, 

TESTl - Cantilever Beam With Axial Load 


A cantilever beam is loaded with a single static compressive loading incre- 
ment. The beam (Figure 3.2-2) is made up of one member, with all degrees of 
freedom constrained at one end and all but two constrained at the end where 
the load is applied. The simple material model is used, and the loading causes 
only elastic displacements. The lowest frequency and buckling factor are ob- 
tained. The displacements, strains and stresses are found to be: 

ui s P/K = -10-^ 

El = ui/L = -10-5 


®1 = Eel = -100 
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Figure 3.2-2 Schematic of TESTl Beam 


The resulting lowest frequency and buckling factor are: 


^lowest = S- = 22.5 

CTT 

KUp 

^cr = 

Agreement between these computed values and independent closed-form solutions 
is exact. 

TEST2 - Simply Supported, Centrally-Loaded Square Plate 


A quarter of the square plate is modeled using symmetry boundary conditions 
(Figure 3.2-3). Four outside beams and four interior diagonal beams are used, 
with dimensions of the beams chosen so as to reproduce the stiffness and mass 
of the plate. One static loading increment is used with the simple constitu- 
tive model in the elastic range. The nonconstrained degrees of freedom are 
shown. 

The theoretical central displacement is: 

U2 = .01160 Pa2/D 
U2 = -3.2428 X 10"6 

The result from the MOMM computer run is: 

U2 = -3.4712 X 10-6 
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a = 8 

D = 9.1575 X 105 


Figure 3.2-3 Square Plate Approximation Centrally Loaded 


TEST3 - Beam With Axial Enforced Displacement (Static) 


A static analysis (Figure 3.2-4) is performed using Walker's viscoplastic ma- 
terial model with twelve loading increments. The properties of Hastelloy X at 
a temperature of 87l“C (1600*F) are used, and the tip displacement is enforced 
at a strain rate of 3.9 x 10-3 sec-1, jhe computer program reproduces the 
experimental results. A plot of the stress-strain curve obtained from the out- 
put is shown in Figure 3,2-5. 



ui 




Au = .0015 
At = .07692 



Figure 3.2-4 Schematic of TESTS 
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Figure 3,2-5 Stress-Strain Response for TESTS 
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TEST4 and TESTS - Beam With Axial Enforced Displacements (Transient) 

Both test cases contain a beam fixed at both ends with a node in the middle of 
the beam (Figure 3.2-6). One end is displaced so that the strain rate equals 
3.9 X 10~3 sec”l. A transient analysis is performed, with TEST4 containing 
Walker's viscoplastic material model and TESTS containing the state-of-the-art 
material model. The viscoplastic materj'al model uses the properties of Hastel- 
loy X at a temperature of 87l“C (1600“F). Figure 3.2-7 shows the displacement 
at the enforced displacement node, as well as the displacement at the center 
node versus time for each model. The results agree exactly with those obtained 
using a simple Euler integration. 



Figure 3.2-6 Schematic of TEST4 and TESTS 



TIME (sec) 


Figure 3.2-7 Displacement History for TEST4 and TESTS 
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3.2,3 List of Input Parameters 


The following list contains the definitions of all variables that are included 
in the input to the Mechanics of Materials Model computer program. 


NGP 

- 

Number of grid points 

GP(I.l) 

- 

X-coordinate of grid point I 

GP{I,2) 

- 

Y-coordinate of grid point I 

GP(I,3) 

- 

Z-coordinate of grid point I 

NB 

- 

Number of beams 

icn.i) 

- 

Grid point at end A of beam I 

IC(I,2) 

- 

Grid point at end B of beam I 

IC(I,3) 

- 

Orientation grid point of beam I 

IC{I,4) 

“ 

Geometry set number of beam I 

IC(I,5) 

- 

Material set number of beam I 

NG 

- 

Number of geometry sets 

BCd.l) 


Width of beam along local y-coordinate 
geometry set I 

BC(I,2) 

“ 

Width of beam along local 2 -coordinate 
geometry set I 

NM 

“ 

Number of material sets 

XMATd.l) 

- 

Young's modulus in material set I 

XMAT(I,2) 

- 

Poisson's ratio in material set I 

XMAT(I,3) 

- 

Coefficient of expansion in material set I 

XMAT(I,4) 

-- 

Zero stress temperature in material set I 

XMAT(I,5) 

- 

Yield stress in material set I 

XMAT(I,6) 

- 

Mass density in material set I 
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TEMP(l) 

SLOPE 


Temperature at initial conditions 
Hardening slope for simple material model 


N - Number of grid points containing degrees of 

freedom 

NI - Number of load increments 

NIP - Number of integration points in each direc- 

tion in each beam 

EE - Static analysis: convergence parameter for 

energy change between two consecutive it- 
erations 

0 Transient Analysis: convergence parameter 
for change in internal energy for the 
adaptive time step calculation 

JJ - Boundary conditions for each degree of freedom 

JJ = 1: nonconstrained 

JJ = 0: constrained 

ICM - Constitutive model indicator 

ICM = 0: simplified material model 

ICM = -1: Walker's elastic-plastic-creep ma- 

terial model 

ICM = 1: state-of-the-art material model 

ITRAN - Transient problem indicator 

ITRAN = 0: no transient analysis 

ITRAN = 1: transient analysis, forces input 

ITRAN = 2: transient analysis, enforced dis- 

placements input 

ISIC - Indicates modulus slope used for frequency 

calculation 

ISIC = 0: initial slope 

ISIC = 1: current slope (not available in 

transient problem) 

ILA - Indicates choice of solving for lowest or all 

frequencies 

ILA = 0: lowest frequency 

ILA =1: all frequencies 
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IBUCK 


ISICB 


DF 

DTEMP(l) 

DTEMP(2) 


Buckling problem indicator 

IBUCK = 0: no buckling analysis 

IBUCK = 1: bucklibg analysis (not available 
in transient problem) 

Indicates modulus slope used for buckling 
calculation 

ISICB = 0: initial slope 

ISICB = 1: current slope 

Applied load or enforced displacement for 
each degree of freedom for an increment 

Temperature increment 

Time increment 
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3.2.4 List of Symbols 


List of Symbols 
Referenced Within Section 3.2 


Symbol 

Description 

Page 

Ml 

Acceleration vector 

3.2-4 


Incremental force vector 

3.2-4 

{vl 

Velocity vector 

3,2-4 

[M] 

Mass matrix 

3.2-4 

DT 

Time step 

3.2-4 

X" 

Eigenvalue 

3,2-4 

[I] 

Identity matrix 

3.2-4 

]xl 

Eigenvector 

3.2-4 

f 

Frequency 

3.2-6 

^cr 

Critical buckling factor 
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3.3 SPECIAL FINITE ELEMENT MODEL 


3.3.1 Literature Survey 

Numerical technology available for use in nonlinear analyses of turbine engine 
hot section components was investigated. An extensive literature survey was 
conducted for MARC by Professor T. J. R. Hughes, Stanford University, as a 
subcontractor, containing 352 references (Reference 1). 

The first topic considered involved recent developments in global solution 
techniques for nonlinear finite element equations, including time discreti- 
zation methods and strategies for nonlinear quasi-static analyses. Literature 
on linear equation solvers was also covered. 

As the second item, the finite element basis and approaches for constructing 
element stiffness matrices were studied with particular emphasis on continuum 
and shell elements. Special elements were also discussed, including elements 
for fracture mechanics applications. 

The numerical treatment of constitutive models, in particular problems asso- 
ciated with computational plasticity, was surveyed, including the effects of 
large strains and rotations. 

Topics such as the self-adaptive mesh refinement associated with error esti- 
mates and error indications, and developments in hardware configurations in 
association with coding strategies were included in the survey. 

A few papers concerning the finite element modeling of hot section components 
were uncovered and reviewed. These papers were mainly concerned with linear 
systems with simplified geometry. 

In the final section of the survey, new developments in the numerical treat- 
ment of contact and friction conditions were studied by virtue of modern math- 
ematical concepts of variational inequalities. 

Overall, the current literature indicates that extensive research and devel- 
opment needs to be carried out on new finite element code concepts in order to 
obtain the significant gains in computational efficiency needed for three-di- 
mensional inelastic analysis of hot section components. 

Solution strategies in nonlinear finite element processes have been given much 
attention in recent years. In particular, a class of computationally efficient 
iterative solution schemes such as the quasi -Newton type techniques has been 
developed for solving nonlinear finite element equations. To avoid reassembly 
and refactorization of the tangent matrix and yet to achieve the quadratic 
convergence properties of the full Newton-Raphson method, line-search, sub- 
space search and secant search techniques have been introduced in the litera- 
ture and have proven useful for certain classes of problems. Such techniques 
are also combined with an automatic load increment size control strategy, usu- 
ally referred to as arc-length type methods. A sophisticated iterative solu- 
tion scheme with an adaptive nonlinear incrementation procedure is one of the 
key ingredients for the solution of the present problem. 
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Possible exploitation of dynamic relaxation and iterative solvers for linear 
equations is suggested in the literature survey. The results reported to date 
show that these methods are perhaps useful for well -conditioned problems, but 
no convincing numerical results are yet reported for ill-conditioned problems 
of engineering importance such as plates and shells, meshes with distorted 
elements, and incompressible problems. Iterative procedures for linear equa- 
tion systems resulting from nondisplacement methods need to be investigated 
from a slightly different point of view in order to obtain some computational 
advantage. It is anticipated that a straightforward implementation of algo- 
rithms available in the literature is not robust enough to meet the present 
purposes. 

In the context of nonlinear dynamic analysis of hot section components, a num- 
ber of time integrators were considered. It was strongly suggested that parti- 
tion and operator splitting methods needed to be investigated in conjunction 
with appropriate automatic schemes to determine the time step size. From such 
a point of view, exploitation of a class of single step algorithms is possible 
and perhaps most appropriate. All the important algorithms that appeared in 
recent publications were covered in the literature survey. Progress subsequent 
to the survey has also been reviewed at MARC. 

It should be noted that no systematic investigation to date has been reported 
on the methodology of adaptive time stepping algorithms. This is in contrast 
to the numerous reports and scientific papers that have been found for auto- 
matic load incrementation in the context of quasi-static finite element analy- 
sis. It seems to be a major task to develop a dynamic transient solution algo- 
rithm capable of handling nonlinear problems in a stable manner and control- 
ling the time increment size adaptively. Only a few references dealing with a 
rather primitive version of such a numerical solution were found in the liter- 
ature survey. 

A number of mathematical contributions noted in the survey were associated 
with a posteriori error estimates and algorithms for adaptive mesh refinement 
based on these estimates. These techniques were investigated only in the 
framework of two-dimensional linear elasticity. 

As summarized above, no existing method was directly applicable to three-di- 
mensional inelastic analyses of turbine engine hot section components. Exist- 
ing numerical technologies, together with the development of new methodolo- 
gies, would be essential ingredients for an efficient finite element proce- 
dure. Due to the three-dimensional nature of the problems in the present 
project, it was anticipated that some sort of iterative approach could sub- 
stantially improve the efficiency of the computational procedures. 

From a computational point of view, algorithmic aspects and coding strategies 
form the most important ingredients for successful numerical simulation of the 
present problem. Indeed, a number of papers addressing this aspect of finite 
element technology v/ere included in the summary. However, the wide variety of 
hardware configurations currently available does not allow us to establish a 
single coding strategy to exploit all of the computing power available. It was 
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observed that a basically sound and transparent code could be modified to im- 
prove its performance in a specific hardware environment, but if a code is de- 
signed to maximize its performance in a specific hardware and software config- 
uration, then it would become extremely difficult to effectively transplant 
such a code to different systems. 

3.3.2 Formulation Development 

3.3.2. 1 Introduction 

A finite element solution strategy was designed with particular emphasis on 
three-dimensional inelastic analyses of turbine engine hot section components. 
Key ingredients employed in this strategy are mixed variational formulations 
and their iterative implementations; linear isoparametric finite element in- 
terpolations with sophisticated integration techniques; advanced techniques in 
computational plasticity, in particular the integration of rate independent 
constitutive equations; and a class of single step, second order time integra- 
tors. 

For the spatial discretization, a version of the Hellinger-Reissner variation- 
al formulation for solid mechanics is utilized as the basic variational state- 
ment of the problem, where the displacements and strain components are taken 
as the field variables. A linear Lagrangian finite element basis is used for 
the interpolation of these variables as well as for the stress-strain law. The 
radial return concept plays a central role in the computational plasticity. 
The incremental iterative solution algorithm is, however, cast in the frame- 
work of mixed finite elements, which results in a different algebraic system 
of equations and hence, different convergence properties, which are generally 
better than those of traditional displacement finite element methods. 

The transient algorithms are looked at in a weighted residual manner, i.e., 
the time-space field is split in a logical fashion with the implicit and un- 
conditionally stable nature of embedded time finite element discretization 
being maintained, but the evolutional nature of the original problem pre- 
served. One of the major exercises in the formulation development is to con- 
struct a reasonable engineering criteria to determine the optimal step size at 
each time increment. 

Thus, for the class of problems stated in the definition oF the task, the spa- 
tial discretization method and the iterative procedure for linear and nonlin- 
ear problems have been firmly established. The current procedures still leave 
open certain possibilities for further improvement of convergence properties, 
if higher order solution schemes such as the conjugate gradient method are 
utilized. The transient algorithm, as well as its theoretical background in 
conjunction with present spatial discretization techniques, needs further in- 
vestigation in order to fully utilize the advantages of the present mixed it- 
erative techniques in the context of nonlinear dynamics. 
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3. 3. 2. 2 Mixed Forms and Iterative Solutions 

It is convenient to discuss the framework of the finite element methodology 
used here in the rather general context of a three-dimensional continuous 
body. The deformation is assumed to be small. For clarity, the development is 
presented in terms of classical elastic conditions. The generalization to in- 
cremental inelastic analysis is straightforward and involves the use of an 
appropriate constitutive relationship together with incremental forms of the 
variational equations. 


Consider a deformable body n in three-dimensional physical space, of which the 
boundary afi is sufficiently smooth. Motion and deformation of the body is as- 
sumed small. The deformation and stress history of the body is characterized 
by three field variables: the displacement u, the strain e, and the stress a. 
Using lower case subscripts to denote rectangular Cartesian components of vec- 
tors and tensors with respect to a fixed spatial reference frame, the govern- 
ing differential equations are 


and 




ij,i " 

(3.3-1) 

°ij " “^idkl ^kjf 

(3.3-2) 

= 1/2 (u. . + u. .) 

1 J } 1 

(3.3-3) 


where p is the density of material, a is the acceleration of the body given as 
a time derivative of the displacements, and D is a fourth order tensor which 
describes the material response at a given stress and strain state. The vector 
f is the loading function due to the body force. 


For a given initial state, a set of boundary conditions; 
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completes the classical statement of the problem. Following the derivation 
given by Zienkiewicz and Nakazawa (Reference 2), the first variations of the 
Hu-Washizu variational principle are obtained via the Galerkin method of 
weighted residuals: 
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and 
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dv = 1/2 





(3.3-8) 


where * denotes a virtual quantity. In the above equation the virtual dis- 
placement is assumed to satisfy the homogeneous displacement boundary condi- 
ti on . 


The above statement can be used as a basis for construction of a finite ele- 
ment procedure. One of the major advantages of this form over the conventional 
displacement approach is the explicit presence of stress and strain in the 
variational form and thus in the finite element equations. 

This is the main theme of the present fornwlation development and will be dis- 
cussed in the following sections. 


Use will also be made of the classical virtual work 


Note here that equations (3.3-2) and (3.3-3) are 
variational statement. 


Direct discretization of the Hu-Washizu principle 
equation of the form: 
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embedded implicitly in the 
leads to a finite element 
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with u, e and a being the vector of nodal variables associated with the dis- 
placement, the strain and the stress respectively. 


Eliminating the stress terms algebraically, a finite element form similar to 
the Hellinger-Reissner principle is obtained as: 
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+bq-1d 
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from which an iterative procedure is constructed. 
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Using the discretized form of virtual work statement, equation (3.3-9), the 
above form is modified. 


K +BQ-1d 

+dT(Q-1)TbT _d 

with K being the usual finite element stiffness matrix. The recursive form 
used here is: 

Ku"^^ » F + Ku" - BQ’^De" (3.3-13) 

= (Q'^)V u"''^ (3.3-14) 

which is in principle identical to an iterative method proposed by Loubignac 
(References 3 and 4) and investigated further by Cook (Reference 5). It should 
be noted that this procedure is used extensively in optimization theory (Ref- 
erence 6) for a class of minimization problems with more explicit equality and 
inequality constraints. In finite element computations, solution of incom- 
pressible problems has been attempted using this class of iterative tech- 
niques, which results in an algorithm of Uzawa (Reference 7). A recent pub- 
lication by Fortin and Glowinski (Reference 8) covers the theory of the it- 
erative solution for a wide range of constrained problems in mechanics. 
Nakazawa, et al (Reference 9) show that the high order methods discussed in 
Reference 8 improve the convergence properties of this solution strategy sig- 
nificantly. The first attempt to unify the concept of an iterative solution 
for mixed finite element methods is reported in Reference 10 where another al- 
ternative economical way Tor solving mixed finite element equations is dis- 
cussed. 

In the computational procedure, the algebraic form is treated in an incre- 
mental way, i.e., 

u"^^ = u" + K"^ (F - BQ~^Ds") (3.3-15) 

and 

= (Q'^)^B^u”'’^ (3.3-16) 

The stress is recovered directly from the strain, and therefore varies in a 
similar way as the strain, i.e.. 


;f ^ Ku( 
: 0 ) 


(3.3-12) 


(3.3-17) 

This completes the discussion of discretization and solution procedures for a 
general class of problems in solid mechanics. 


3.3-6 



In the present algorithm, the solution is initialized by the conventional dis- 
placement stiffness equations and this stiffness array is used as the pivot in 
the subsequent iterative solution. Because all the necessary conditions for 
existence, uniqueness and stability are satisfied, the quality of the con- 
verged solution is dictated by the strain interpolation. The displacement 
plays a rather insignificant role other than as a preconditioner to the iter- 
ative approach. This has been noticed in the early papers by Cantin, Loubignac 
and Touzot (Reference 4). 

Using appropriate diagonal ization techniques for the matrix Q, the recovery of 
nodal strain components does not involve a matrix inversion operation, and no 
significant additional effort is required to compute and iterate on this quan- 
tity. 

Once the matrix Q is diagonalized, the procedure to recover the nodal stress 
becomes extremely simple requiring the evaluation of the constitutive law at 
nodes where the strain is calculated. 

The iterative solution procedure is readily applicable to a class of nonlinear 
material problems such as rate-independent plasticity and involves evaluating 
the constitutive law at each iteration cycle. This approach results in a 
scheme similar to the Newton-Raphson method as used in finite element dis- 
placement analyses. The only difference occurs in the procedures to evaluate 
the residuals. Compared with the conventional displacement method, the number 
of operations required for formation of the residual vector at each iteration 
is reduced somewhat because the number of nodal points is usually far less 
than the number of integration points in a given finite element mesh. This 
saves time in the constitutive calculation. 

The iterative procedure developed here provides in principle a powerful ve- 
hicle to study large scale inelastic analysis problems in three dimensions. 
The largest array appearing in this calculation is the same as that in the 
conventional displacement method; however, the strain and the stress here are 
evaluated at the nodes, and hence a better approximation for these field vari- 
ables is obtained. This property, in turn, reduces the necessity of exces- 
sively refining the finite element mesh to obtain accurate stress and strain 
fields. 

Numerical instability is often encountered in stress and strain fields calcu- 
lated from a finite element displacement solution; this is often seen as os- 
cillatory behavior in the numerical approximations and may lead to inaccurate 
inelastic response of the discretized model. The present approach, being sta- 
ble in terms of both displacement and strain, eliminates the possible occur- 
rence of such numerical problems. 
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3. 3. 2. 3 Element Technology 


A family of linear Lagrangian finite elements has been used in accomplishing 
the present task (Reference 11), but the element library is independent of the 
global solution strategy and higher order elements can be added. 

In two-dimensional and axisymmetric cases, four-noded quadrilateral elements 
with nodal points located at corners are used, with the displacement and 
strain components taken as the primary variables. The constitutive equation 
(3.3-2) is evaluated at the same nodes, and parameters associated with it are 
stored together with the stress components. Eight-noded isoparametric brick 
elements are utilized in three dimensions, and the same mechanism of repre- 
senting all the quantities at the nodes is used. For the analysis of plates 
and shells, a four-noded Lagrangian element is used incorporating the trans- 
verse shear terms into a modified variational formulation of Reissner-Mindl in 
type. 

In order to implement the analysis procedures established in the previous sec- 
tion, a number of terms must be integrated to form coefficient matrices. As 
discussed in the context of penalty finite elements for incompressible prob- 
lems, due care must be exercised in order to obtain a stable finite element 
approximation. The integration options for linear quadrilateral elements are 
four-point and single-point Gaussian quadratures, as well as the four-point 
trapezoidal rule, as shown in Figure 3.3-1. For the three-dimensional solid 
elements, eight-point integration rules that correspond to the two-dimensional 
four-point schemes are employed. 

Evaluation of the stiffness matrix K associated with the displacement formu- 
lation for the two- and three-dimensional elements may be done using either a 
standard or a selective reduced integration procedure. In the selective inte- 
gration procedure, direct (or normal) strain components are evaluated at four/ 
eight Gaussian quadrature points, whereas the shear components are dealt with 
at the centroid of the element. This simple "trick" greatly improves the be- 
havior of the elements, particularly in bending. In order to account for the 
isoparametric distortion of elements, the terms in the stiffness equations are 
represented with respect to a local element Cartesian coordinate system, the 
definition of which is given by Nagtegaal and Slater (Reference 12) and shown 
in Figure 3.3-2. 

The load vector (including the surface traction term) is integrated using full 
Gaussian quadrature. The residual vector which appears in Equation (3.3-15) is 
also fully integrated. 

In the strain recovery phase, the use of reduced (single-point Gaussian) quad- 
rature has been found stable and accurate, but the quality of the displacement 
solution often deteriorates after a few iterations for a given increment. The 
use of the trapezoidal rule which naturally results in the lumped (diagonal) 
form of the Q matrix is found to be optimal. For certain cases, the quality of 
the displacement solution is preserved when this integration procedure is 
used, without any loss of accuracy in stress and strain approximations. 
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Figure 3.3-2 Coordinate Transformation for Two-Dimensional Elements With Iso- 
parametric Distortion 



For the shell element, an adjustable form of the selectively reduced integra- 
tion technique has been applied to the integration of displacement stiffness 
terms, to avoid singularity and numerical locking. The same options are made 
available for the strain recovery at nodal points. To recover the strain at 
nodes, a local Cartesian coordinate system common to all the elements joined 
to the node is required. Vector and tensor components, initially defined with 
respect to the shell element coordinates, are transformed to the local Car- 
tesian system using appropriate operations. 

3. 3. 2. 4 Inelastic Constitutive Models 

Three inelastic constitutive models have been implemented into the MHOST pro- 
gram. In order of increasing sophistication, these models are as follows: 

Approximate Constitutive Model 


This model approximates the global stress-strain behavior with a nonlinear 
pseudo-elastic model. The model is described in detail in a report by 
Cassenti (Reference 13). A secant modulus procedure is used for the gen- 
eration of stress increments and constitutive equations. This simple model 
does not allow for numerous important effects actually occurring in the 
cyclic loading of hot section components. Nevertheless, it definitely has 
some usefulness in obtaining some indication of the degree of nonlinearity 
in simple load simulations. 

State-of-the-Art Constitutive Model 


This model includes the approaches most commonly used for inelastic analy- 
sis of hot section components. It combines state-of-the-art time-indepen- 
dent plasticity with isotropic (and in the near future, also kinematic and 
combined) strain-hardening, a classical creep model in series with the 
plasticity model and a thermal expansion capability. In addition, elastic 
and inelastic properties can be made temperature-dependent. The material 
may contain initial anisotropic behavior for the elastic properties, yield 
surface definition and coefficient of thermal expansion. 

In the plasticity model, the radial return method is used for stress re- 
covery and the tangent modulus approach is used to generate the stress- 
strain law. In the creep model, a simple explicit procedure is applied for 
time integration of the creep strains. Both methods can also be regarded 
as state-of-the-art methods for implementation of such models. 

The state-of-the-art model includes most of the important effects occur- 
ring in cyclic loading of hot section components. However, some weaknesses 
exist in load reversal situations, in particular with respect to interac- 
tion between plasticity and creep effects; if these effects are very im- 
portant, the advanced model should be employed. 
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Advanced Constitutive Model 


The integrated creep-plasticity model orginally developed by Walker (Ref- 
erence 14) and later modified by Cassenti (Reference 13) was implemented 
into the HOST program. The model is described in detail in the report by 
Cassenti (Reference 13). Although both time- independent plasticity and 
creep effects can be distinguished in the model, definite interactions 
have been included on the single-point level. In addition, reversed load- 
ing effects are included more accurately in this model than in the state- 
of-the-art model . 

An explicit subincrement procedure with optimal time step selection is 
used in the stress recovery procedure for this model. The elastic (but 
temperature dependent) stress-strain law is used for the formulation of 
the stiffness matrix; hence, the iterative algorithm for this model is, in 
fact, the modified Newton method. 

The advanced constitutive model allows accurate material modeling for hot 
section components. It should be noted, however, that the complexity/ 
generality of the model is also its major drawback. Almost a dozen tem- 
perature dependent material constants need to be determined experimentally 
before the model can be applied. 

3. 3. 2. 5 Time Integration 

For a class of nonlinear problems to be dealt with in the present task, dynam- 
ic effects were taken into account. A recursive form is derived and implement- 
ed here, which is conceptually based upon the embedded time (temporal) finite 
element concept. 

The solution of an evolution problem lies in a four-dimensional time-space. 
Therefore, in the formulation it is necessary to take temporal as well as spa- 
tial variations into account using a four-dimensional basis for the discreti- 
zation. However, as demonstrated in Section 3. 3. 3. 2, the differential opera- 
tors governing the evolution of nonlinear structural mechanics can be split 
apart and the dependency of the solution upon each independent variable can 
now be treated separately. Such a separation of variables is of some impor- 
tance in dynamic situations, because the differential equation consists of 
different characteristics with respect to each independent variable, that is 
hyperbolic in time and elliptic in all spatial dimensions. 

To preserve the evolutional nature of the original problem and to take into 
account the hyperbolic nature of temporal operator, the discretized system of 
the embedded time finite elements is used recursively from one time level to 
the next.^ Hence, the equivalent form is derived by integrating the spatially 
discretized form of finite element equations in time. It is noted that a non- 
symmetric weighting function is used in this procedure due to the well-known 
conditional stability of finite element procedures for hyperbolic systems. 
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Using the two-noded time element and an appropriate nonsymmetric weighting 
function, the algebraic system becomes identical to the generalized Newmark-s 
finite difference expressions. 

In each time element, the displacement increment obtained by the integration 
scheme is fed into the same discrete strain-displacement equation (3.3-14) to 
obtain the first correction vector, and then the residual vector is updated to 
produce the new displacement increment. The residual at the end of each time 
increment is carried to the next time level so that the overall equilibrium is 
satisfied globally in the sense of embedded time finite element methods. 

The size of the time element is controlled adaptively with respect to an ener- 
gy weighted average of the relative phase errors, based upon the contributions 
of strain and kinematic energy in the deformation represented by a time ele- 
ment. 

3, 3. 2. 6 Eigenvalue Extraction 

Eigenvalue extraction is utilized in the MHOST program to obtain both the dy- 
namic model frequencies and collapse load estimates. In addition to the eigen- 
values, the eigenvectors are also obtained. 

The subspace iteration technique has been implemented in the MHOST program. 
Using this method, the global stiffness and mass matrices are transformed into 
a subspace. A threshold Jacobi method has been used to obtain all of the 
eigenvalues in the subspace. This method has been found to converge very 
quickly. The program iterates forming new subspaces until convergence is ob- 
tained. 

This method has been developed for the extraction of a large number of modes 
in a large system. This is typically what is required for modal dynamic use. 
The method has not been applied as frequently to the calculation of buckling 
modes, where only a few modes are required. It is anticipated that this is an 
area where future work may be fruitful. 

An important issue in the successful convergence of this method is in the 
choice of the initial trial vectors. Currently, a simple scheme is being used, 
where the degrees of freedom with smallest K/M ratio are given the greatest 
weight. 

Currently the lumped mass matrix is being utilized. 

3.3.3 Program Development 
3. 3. 3.1 Introduction 

A finite element code (to be referred to as the MHOST program) has been de- 
veloped by MARC incorporating the formulation and solution strategies pre- 
viously described. The framework of the code is built on the well established 
foundation of finite element displacement coding for nonlinear structural 
analysis. 


3 . 3-12 



The mixed interpolations and the iterative procedures are built into the MHOST 
program as additional operations to the standard procedure, and no major 
changes are required inside of the finite element computations. One of the ma- 
jor differences, however, from the user's point of view is that all informa- 
tion is available at the nodal points rather than at the element integration 
points. 

The current version of the MHOST program is essentially a vehicle to explore 
the capability of the present strategyaand no attempt has been made to opti- 
mize the computational efficiency at the coding level. The program is written 
in a clear manner in order to achieve high productivity in the formulation and 
program development. 

In this section, a technical note is provided on the overview and control 
structure of the MHOST program, the element library and the nonlinear solution 
capability of the system. 

3. 3. 3. 2 Overview and Control Structure 

The MHOST program is written in FORTRAN IV with commonly accepted options, and 
has been tested on the PRIME Primes 18.2 FTN compiler and on the IBM/CMS 
FORTRAN H extended compiler with the optimization level 2. The advantages of a 
virtual storage system are exploited by storing most of the information in the 
core memory, and no special file input/output operations are utilized. Hence, 
the code is portable to other virtual memory systems, with little conversion 
effort necessary. The code can be run either interactively or in a remote 
batch mode. 

The main functions of the MHOST code are illustrated in Figure 3.3-3, which 
depicts a number of processors devoted to specific operations. In the follow- 
ing paragraphs, the functions of these processors will be briefly discussed. 

A user-friendly, free-format input processor is attached which reads in the 
data from the main input channel. The input data consist of PARAMETER DATA, 
MODEL DATA and INCREMENTAL DATA. The PARAMETER DATA specify the size of the 
model and select the options for the integration and the types of loading. The 
MODEL DATA provide the definition of the finite element model including the 
boundary conditions and load data for the first increments. The INCREMENTAL 
DATA section provides additional loading data and boundary conditions for each 
load step of the incremental solution. 

Data of different types within the three major data blocks are identified 
using keywords. The keywords available for use in the PARAMETER, MODEL, and 
INCREMENTAL data blocks are listed and defined in Section 3.3.5. 
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Figure 3,3-3 The MHOST Program Main Control Flow 
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The input processor passes information to the main finite element processor in 
which the iterative solution for mixed finite element equations is carried 
out. The flow chart of this section of the MHOST code is given in Figure 
3.3-4, which also includes additional operations to the conventional displace- 
ment solution indicated by bold lines. 

During an increment, the initial solution of the iteration for the mixed pro- 
blem is carried out in exactly the same manner as the displacement method, 

with the converged stress and strain field of the previous increment taken as 
the initial state to calculate the tangent array. The mixed interpolation then 
takes over the process, and the tangent array becomes part of the strategy. 
The nodal strain is recovered from the displacement solution, and then the 

constitutive equation is integrated at the nodal point level to obtain the 
nodal stress array. The new residual vector is formed using this new state of 
stress. The iteration is repeated until the norm of the residual vector be- 
comes less than the prescribed tolerance. 

The size of the new incremental load is adaptively calculated when the auto- 
matic incrementation routine is invoked. Otherwise, the input data for the 

next increment are loaded and the above procedure is repeated until the stop 
flag is set. 

The output processor generates a report printout on the line printer channel 
at the end of every increment. When the restart file and the plot file are 
requested, this processor generates these files and catalogues them in the 
system. 

As dictated by the mixed formulation of the problem, all the information pro- 
duced by the MHOST code is associated with nodal points. In the options to 
generate the line printer file, however, the element integration point infor- 
mation is made available using the values interpolated by the shape function 
from the nodal values. 
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Figure 3.3-4 


Detailed Flow Chart for Iterative Solution for Mixed Problem 
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3. 3. 3. 3 Element Library 

The control structure of the code is independent of the element type used in 
the analysis described in the formulation development Section 3. 3. 2. 2. The 
operations associated with the element technology are coded in the element 
library subroutines. This class of subroutines consists of the shape func- 
tions, the derivatives of shape functions with respect to the global coordi- 
nate system and the element Cartesian coordinate transformations. These in- 
dividual components are merged into the macro subroutines which assemble ma- 
trices used in the solution procedure. These macro subroutines are constructed 
for every type of geometric feature of the model, i.e., plane stress, plane 
strain, axisymmetric and, three-dimensional continuum. When the MHOST program 
is executed, the upper level control routines select the appropriate macro 
element library subroutines; and then each macro subroutine, in turn, selects 
the lowest level element library subroutines. 

In the current version of the MHOST program, the bottom level element library 
consists of only two element types, i.e., a four-noded two-dimensional element 
and an eight-noded three-dimensional element. The structure of this code al- 
lows the user to enhance the element library without major difficulties. 

The macro level subroutines consist of the lumped mass matrix routine to form 
Q, and the displacement strain matrix routine to form B. No element level 
operation is necessary for the stress recovery operations. For shells, the 
geometrical features of the problem require some coordinate transformations 
which are defined only at the element level. 

3. 3. 3. 4 Nonlinear Analysis Capabilities 

The following geometric types are included in the MHOST program; 

1. truss in three dimensions, 

2. plane stress, 

3. plane strain, 

4. axisymmetric solid, 

5. three-dimensional solid, and 

6. shell in three-dimensions. 

The nonlinear constitutive equations built in this code are evaluated once at 
each iteration loop when the nodal stress is recovered after the strain-dis- 
placement solution. 

The material laws discussed in the formulation development Section 3. 3. 2. 2 are 
readily built into the present version of the MHOST code leaving the possi- 
bility of major modification by means of user subroutines. 
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3.3.4 Program Val idation/Ven'fi cation 


The validation/verification of analysis capabilities is an important step in 
the computer program development process. Validation/verification of the MHOST 
finite element program is a joint MARC/Pratt & Whitney effort, and involves 
comparing MHOST results with those obtained from closed form solutions, the 
MARC General Purpose Structural Analysis program, and hardware tests. 

The simplest validation cases involved subjecting several element types to 
uniaxial loading. The finite element mesh and applied loads for a plane stress 
block are depicted on Figure 3.3-5. After three increments of loading into the 
plastic range, MHOST and MARC results agreed to three significant figures. The 
mesh shown on Figure 3.3-5 was also analyzed using the four node plane strain 
quadrilateral element and agreement between results from the two finite ele- 
ment codes was also very good (three significant figures). In addition, a 
similar mesh of axi symmetric elements was exercised for the case of uniaxial 
tension. After three increments of loading, results from both codes agreed to 
two significant figures, except for the plastic strains. Here, MHOST pre- 
dictions were about 2 % smaller than MARC values. 

The three-dimensional block shown on Figure 3.3-6 was placed in uniaxial ten- 
sion by prescribing displacements in the x-direction. The MHOST results exhib- 
ited homogeneous uniaxial behavior in the plastic range after three increments 
of prescribed displacement, and numerical values were in excellent agreement 
with a closed form solution. 

A simple linear temperature analysis case was created by constraining the ends 
of the mesh shown on Figure 3.3-5 and subjecting it to a uniform increase in 
temperature. MHOST results were identical to the closed form solution. 

The more complex case of the elastic response of a cantilever beam subjected 
to a moment load was also analyzed using the MHOST program. The rectangular 
finite element mesh and loading data are shown on Figure 3,3-7. Five Loubignac 
iterations were performed and the results obtained using various numerical 
integration procedures are presented in Table 3-1. Run A which employs selec- 
tive stiffness integration and strain-displacement integration by the trape- 
zoidal rule shows the best (essentially exact) results. 

Highly irregular element shapes are often used in models of gas turbine com- 
ponents in order to minimize problem size and, hence, control analysis costs. 
Since the performance of most quadrilateral and hexahedron elements deterior- 
ates with increasing departure from parallelepiped shapes, the shape sensiti- 
vity of these elements must be quantified in order to evaluate the accepta- 
bility of meshes proposed for practical applications. The effects of Loubignac 
iteration on the shape sensitivity of the quadrilateral plane stress element 
were examined using a cantilever beam modeled with skewed meshes and subjected 
to a pure moment load as shown on Figure 3.3-8. Results were obtained for 
three values of the mesh skewing parameter 0(0®, 22,5", 45") defined on the 
figure. 


3.3-18 



Figure 3.3 
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-6 Finite Element Mesh and Prescribed Displacements Used for the 
Validation of Three-Dimensional Elements 
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Figure 3,3-7 Cantilever Beam Problem (Plane Stress) 


Table 3-1 

Cantilever Beam; Loubignac Iteration Comparison 



Exact 

Run A 

Run B 

Run C 

Tip Displacement (in) 

2.4 

2.3988 

3.4951 

2.2699 

Stress (psi) 

6000 

6005.1 

4368.1 

5710.7 


Sti ffness 
Integration 

Run A Selective 

Run B Full Gaussian 

Run C* Full Gaussian 


*The results are identical to the full Gauss stiffness matrix integration in . 
conjunction with the reduced integration strain recovery integration. 
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Figure 3.3-8 Cantilever Beam: Shape Sensitivity Study 


The behavior of the tip displacement and a typical extreme fiber bending 
stress are plotted on Figure 3.3-9 as functions of the number of iterations. 
Values for the cases of full and selective stiffness integration are shown on 
the figure (strain-displacement integration was performed using the trape- 
zoidal rule in both cases). Displacements and stresses approach exact values 
as the number of iterations is increased and the improvement in these results 
is dramatic in the first few iterations. Hence, Loubignac iteration is an ef- 
fective technique for minimizing the deterioration in accuracy associated with 
the use of irregular meshes. The positive effects of using selectively reduced 
element stiffness integration are apparent; two to four fewer iterations are 
required to obtain results equivalent to full element stiffness integration 
values. 

The MHOST four node shell element was used to analyze a simply supported 
square plate under the action of a uniform pressure. The finite element model 
including geometry, loads and boundary conditions is shown on Figure 3.3-10. 
An elastic analysis with five iterations was performed. The thickness/length 
ratio for the plate is very small (0.005) so comparisons of MHOST results with 
the classical thin plate solution in Reference 15 are appropriate. The exact 
and computed values for lateral displacement and My at the center of the 
plate as well as the values for My at the x = 1.0, y = 0.0 location are 
shown in Table 3-1 1. The MHOST results at the center of the plate are excel- 
lent (errors less than 2 percent) and an adequate representation of the zero 
moment at mid-edge is also obtained. 
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Figure 3.3-9 Convergence of Iterative Procedure for Elastic Cantilever Beam 
Modeled With Skewed Meshes 
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Figure 3.3-10 MHOST Model of Simply Supported Plate 


Table 3-II 

Simply Supported Plate: Uniform Pressure Load 


Item 

Exact 

MHOST 

Center Displacement (in) 

-7.093 

-7.231 

My at X = 1.0, y = 0,0 (in-lb/in) 

0.00 

-0.139 

My at Center (in-lb/in) 

-1.916 

-1.932 
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The MHOST program has been used to determine the elastic-plastic response of a 
thick cylinder under the action of a uniform internal pressure p-j. Figure 

3.3- 11. The cylinder has internal radius r=a=1.0 and external radius r=b=2.0, 
with an axial (i.e., z-direction) thickness of 1.0. The cylinder is modeled by 
twenty plane strain axisymmetric elements, each with radial thickness 0.05, 
Figure 3.3-11. Plastic behavior starts at the internal radius and propagates 
radially outward as p-j increases. The material is elastic-perfectly plastic 
with yield stress aY=45000. Other material properties are Young's modulus 
E=30. X 10® and Poisson's ratio v = 0.3. 

Results for this case are shown in Figure 3.3-12, where the axial stress aj 
(normalized) in the elastic region is plotted as a function of the internal 
pressure loading p^ (also normalized). Here, k = Gy/Vs. For completely 
elastic behavior, the axial stress is constant in the radial direction and 
varies linearly with p^- along the straight line in Figure 3.3-12. For elas- 
tic-plastic behavior, the axial stress az is constant only for the elastic 
portion of the cylinder, but varies considerably in the surrounded elastic- 
plastic region. The value of az plotted here is that associated with the 
constant elastic region result. The MHOST program values clearly show very 
good agreement with the theoretical results obtained by Prager and Hodge 
(Reference 16). 

In addition, a plate with a central hole was analyzed via MHOST and results 
were compared with MARC values for identical geometry and loading. The finite 
element mesh for this case is pictured in Figure 3,3-13, where both node and 
element numbering are indicated. The hole of radius 1.0-inch lies in a 6-inch 
by 6- inch square region, with symmetry conditions existing along both the hor- 
izontal and vertical axes, so that only one-quarter of the region need be con- 
sidered. Appropriate displacement symmetry boundary conditions are imposed 
along the left side (i.e., at nodes 28, 30, 25, 20, 22, 24, 35) and along the 
bottom side (i.e,, nodes 15, 18, 10, 3, 6, 9, 31) boundaries. Uniform pressure 
loading is applied along the top boundary, i.e., along side 33-34 for element 
23 and along side 34-35 for element 24. The elastic material properties used 
are as follows: E = 30. x 10®, v = 0,3. The yield stress is 30. x 10^, 
with a piecewise linear stress-plastic strain curve applying after yielding. 

Comparisons of the results are presented in the next two figures. Figure 

3.3- 14 shows the normalized ay stress at the most critical internal radius 

location (i.e., node point 15) as a function of normalized cr|_0Ap external 
load. In addition, corresponding values are also presented for the integration 
point in element 8 lying closest to node point 15. The onom stress is 10 x 
10^ and represents the theoretical external pressure load at which plastic- 
ity is theoretically incipient at node point 15. The node point values associ- 
ated with MARC were obtained from a simple isoparametric bilinear fit of the 
integration point values of element 8. The straight line unconnected to any 
data points represents the theoretical purely elastic response at the inner 

radius. As can be seen from this figure, the ay values agree reasonably well 
between MHOST and MARC; in fact, overall agreement between the two sets of re- 
sults improves as the external pressure loading increases. In a similar man- 
ner, Figure 3,3-15 shows the plastic strain component, -ey*-, varying with 

normalized external load at the same node and integration point locations. 
Again, agreement between MHOST and MARC is quite good, even as the external 
pressure load increases. 
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Figure 3.3-11 Thick Cylinder Geometry Mesh 
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Figure 3.3-14 Normalized ay Stress vs. Normalized External Pressure Loading 


3.3-28 


o 

A 

□ 

o 


-EyPL 

( 10 - 3 ) 

2.0 — 


1.0 I— 


0 

0.8 1 . 


Figure 3.3-15 Plas 
Load 




Several structural problems for which test data is available have been select- 
ed for analysis using the MHOST program. The first of these configurations is 
the NASA benchmark notch fatigue test specimen (Reference 17). 

Elastic analyses of the flat, double-notched uniaxial ly loaded specimen were 
performed using the fine and coarse meshes shown on Figure 3.3-16. In this 
case, the coarse mesh was defined with one-eighth the number of elements em- 
ployed in the fine mesh. The elastic strain results from analyses using the 
two meshes are plotted together with test results in Figure 3.3-17. The micro- 
strain at the surface of the notch for the fine mesh increased from 1705 to 
1818 when three Loubignac iterations were employed. The latter value compares 
favorably with results reported in Reference 17 (1800-1900). Agreement between 
fine mesh values and test results at other locations is also very good. The 
microstrain at the surface of the notch produced by the coarse mesh with three 
iterations (1813) is nearly the same as the iterated fine mesh value and is 
superior to the standard (no iteration) fine mesh value. This result demon- 
strates the effectiveness of the iteration procedure in conjunction with 
coarse meshes and is especially significant because the coarse mesh with it- 
erations used only one- tenth the computer time required to analyze the fine 
mesh without iterations. 

Elastic-plastic analyses of the benchmark notch specimen are now in progress. 
The advanced constitutive model is being employed in these analyses. 

Another MHOST verification problem is related to a combustor application. Con- 
ventional full hoop combustors are subjected to large thermal gradients in the 
radial direction, as well as local streaking in the hoop direction. It is the 
interaction of these effects that causes failure in the form of low cycle fa- 
tigue (LCF) cracking and local oxidation. With these failure mechanisms in 
mind, advanced combustors have eliminated the radial thermal fight by segmen- 
ting the hot wall. This leaves only the local streak temperature variation as 
a driving force for possible LCF failures. The development and optimization of 
segmented configurations requires an understanding of the local material be- 
havior in streak locations. 

Several years ago, the Material Development Group at Pratt & Whitney institut- 
ed an LCF test to rank current and proposed combustor materials. This test 
thermally cycled a flat disk specimen by periodically imposing a small local 
hot spot at the specimen's center. The nature of the test, with a small 
through thickness thermal gradient, causes a bulging or blistering of the hot 
spot toward the flame which is analogous to combustor streak behavior. Analyt- 
ic evaluation (MARC) of both this test and combustor components under streak 
loading confirmed this analogy. In each case, the material in the hot spot 
showed identical deformations. Furthermore, each showed identical amounts of 
damage in the form of creep and plastic strains. It has been concluded that 
the hot spot blister test is a valid simulation of combustor streak behavior. 
Figure 3,3-18 shows a MHOST finite element solid model of a 1/4 segment of the 
flat circular specimen. The model consists of 320 solid elements and 615 
nodes. Validation/verification analyses are currently in process. 
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Figure 3.3-16 Three-Dimensional NASA Benchmark Notch Specimen 


1.9 

1.13: 

1.7l> 

1 . 6 - 
1.5- 
1.4- 
1.3- 
1 . 2 - 
1.H 
1 . 0 - 
0.9- 
0.8 -I 

0.7 


□ 

■f 


• STRAIN GAGE DATA (NASA CR-1 65571 ) 
O FINE MESH (NO ITERATIONS) 

f FINE MESH (3 ITERATIONS) 

0 COARSE MESH (3 ITERATIONS) 








— ( 


1.0 2.0 3.0 4.0 5.0 6.0 

DISTANCE FROM NOTCH (mm) 


7.0 


8.0 


9.0 


Figure 3.3-17 Comparison of Calculated and Measured Elastic Strains for the 
NASA Benchmark Notch Specimen 
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Figure 3.3-18 Simulated Combustor Liner Specimen (One-Quarter Section) 


A sophisticated MHOST finite element model of a low pressure turbine vane has 
been developed to demonstrate the program's applicability to this type of pro- 
blem. The initial model shown on Figure 3.3-19 contains 424 four-noded shell 
elements which is consistent with models that would normally be employed in 
NASTRAN or MARC analyses. However, the enhanced capability of the MHOST code 
(Loubignac type iteration) should allow the ultimate use of a much coarser 
model. The loading environment consists of both thermal and aerodynamic pres- 
sure loads. Since actual vanes of this type are clustered in groups of three, 
cyclicly symmetric boundary conditions are used at the inner vane platform to 
simulate the presence of the additional vanes. Elastic stress/strain distribu- 
tions predicted by the MHOST analysis will be compared directly to strain gage 
measurements taken on actual vanes. In addition, a creep analysis will be per- 
formed to simulate the 69 hours of hot time accumulated at takeoff power dur- 
ing an in-house endurance test of this vane. The analytical creep predictions 
will be compared to deflection data from the actual test. 



Figure 3.3-19 Turbine Vane Finite Element Model 


3.3.5 Input Data Structure for the MHOST Program 

The input data for the MHOST program is divided into three major blocks: 
PARAMETER DATA, MODEL DATA, and INCREMENTAL DATA. Data of different types 
within the data blocks are identified using keywords as discussed in Section 
3. 3. 3. 2. The keywords are listed and defined in this section to provide the 
potential user of MHOST with an overview of the structure and content of pro- 
gram input. 

The PARAMETER DATA serve to specify dynamic array allocations and to control 
the MHOST module execution sequence. The keywords associated with this data 
block are listed below: 


KEYWORD 

SIGNAL TO MHOST PROGRAM 


*ANIS0TR0PY 

Anisotropic material data will be input 


*B0LINDARY 

Displacement boundary conditions will be input 


*BUCKLE 

Perform a buckling analysis 


*CONSTITUTIVE 

Use constitutive model (simplified, conventional, 
or advanced) indicated by parameter value 

*CREEP 

Perform a creep analysis 


*DISTRIBUTEDLOAD 

Distributed loads will be input 


*DUPLICATENODE 

Duplicate nodes will be input 


*DYNAMIC 

Perform a transient analysis 


*ELEMENTS 

Use element type indicated by parameter value 


*F0RCES 

Nodal forces will be input 


*LOUBIGNAC 

Use Loubignac iteration options indicated by 
rameter values 

pa- 

*M0DAL 

Perform free vibrations analysis 


*N0DES 

Upper bound to number of nodes indicated by 
rameter value 

pa- 

*0PTIMIZE 

Optimize band-width of stiffness matrix 


*PERI0DICL0ADING 

Periodic loads will be input 


*PRINTSETS 

Upper bound to number of print sets given by 
rameter value 

pa- 
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KEYWORD 

*REPORT 

*RESTART 

*SCHEME 

*STRESS 

*TANGENT 

*TEMPERATURE 

*THERMAL 

*TRANSFORMATIONS 

*TYING 

*UCOEFF 

*UDERIV 

*UFORCE 

*UHOOK 

*UPRESS 

*UTEMP 

*UTHERM 


SIGNAL TO MHOST PROGRAM 

Print interval given by parameter value 

A previous analysis is to be continued 

Time integration operator identified by parameter 
values 

Upper bound to number of stress boundary condi- 
tions given by parameter value 

Perform modified Newton-Raphson iteration after 
cycle given by parameter value 

Nodal temperatures will be used in calculations 

Temperature dependent material properties will be 
used 

Upper bound to number of coordinate transforma- 
tions given by parameter value 

Number and form of tying equations given by pa- 
rameter values 

User subroutine UCOEFF will be provided 
User subroutine UDERIV will be provided 
User subroutine UFORCE will be provided 
User subroutine UHOOK will be provided 
User subroutine UPRESS will be provided 
User subroutine UTEMP will be provided 
User subroutine UTHERM will be provided 
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The finite element model is defined with the MODEL DATA. For example, the fi- 
nite element topology and nodal coordinates, as well as boundary conditions 
and loads, are specified in this data block. The MODEL DATA keywords are iden- 
tified below: 


KEYWORD 

DEFINITION 

*B0DYF0RCE 

Body force data 

*B0UNDARY 

Prescribed displacement data 

*C00RDINATES 

Nodal coordinate data 

*DISTRIBUTEDLOAD 

Distributed load data 

*DUP LI CATEMODE 

Duplicate node data 

*ELEMENTS 

Element connectivity data 

*END 

End of MODEL DATA 

*F0RCES 

Nodal force data 

*INCREMENTS 

Maximum number of increments 

*ITERATIONS 

Iteration control data 

*P0STFILE 

Postprocessing file control data 

*PRINT0PTI0N 

Output control data 

*PROPERTIES 

Element properties data 

*SAVE 

Restart control data 

*ST0P 

End of analysis 

*STRESS 

Stress boundary condition data 

TEMPERATURES 

Nodal temperature data 

TRANSFORMATIONS 

Coordinate transformation data 

TYING 

Tying equation data 

*W0RKHARD 

Equivalent stress - equivalent plastic strain data 
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The INCREMENTAL DATA are used to specify the loading history, modifications to 
boundary conditions, and initial conditions for transient response calcula- 
tions. Any number of INCREMENTAL DATA blocks may be included in an analysis. 
The following MODEL DATA keywords may be used in an INCREMENTAL DATA block: 


*B0DYF0RCE 

*B0UNDARY 

*DISTRIBUTEDLOAD 

*END 

*F0RCES 

*ITERATIONS 

*P0STFILE 

*PRINT0PTI0N 

*SAVE 

*ST0P 

TEMPERATURES 

TYING 


The additional keywords available for use in the INCREMENTAL DATA blocks are 
defined below: 


KEYWORD 

*ACCELERATION 

*AUTOINCREMENT 

*DISPLACEMENT 

*ENDINITIALCONDITION 

*INITIALCONDITION 

*PERIODICLOADING 

*PR0P0RTI0NAL 

TIME 

TELOCITY 


DEFINITION 


Initial node accelerations data block (transient 
analysis) 

Automatic load increment data 

Initial nodal displacements (transient analysis) 

End of initial condition data (transient analysis) 

Start of initial condition data (transient analy- 
sis) 

Periodic loading data 
Proportional loading data 

Time increment data for creep and transient anal- 
ysis 

Initial nodal velocities (transient analysis) 
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3.3.6 List of Symbols 


List of Symbols 
Referenced Within Section 3.3 


Symbol 

Description 

Page 


Acceleration components 

3.3-4 

B 

Matrix which relates nodal stresses to nodal forces 

3.3-6 

D 

Stress-strain matrix 

3.3-5 

f^ijkje 

Stress-strain tensor 

CO 

• ' 

CO 

1 

E 

Elastic modulus 

3.3-20 


Body force components 

1 

CO 

• • 

CO 

F 

Nodal force vector 

3.3-5 

h 

Thickness of model 

3.3-20 

K 

Conventional displacement method stiffness matrix 

3.3-6 

M 

Moment 

3.3-21 

"j 

Direction cosines of normal to surface 

3.3-4 

Pi 

Internal pressure 

3.3-24 

Q 

Shape function inner product matrix 

3.3-5 

R 

Radius 

3.3-24 

^i 

Surface tractions 

3.3-4 

u 

Displacement vector 

3.3-4 

x,y .2 

Coordinates 

3.3-9 

5 

Displacement component at a point 

3.3-19 

c 

Nodal strain vector 

3.3-5 

^ij 

Strain tensor 

3.3-4 

n»^ 

Isoparametric coordinates 

3.3-9 
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List of Symbols 

Referenced Within Section 3.3 (continued) 


Symbol 

Description 

Page 

e 

Mesh skew angle 

3.3-18 

V 

Poisson's ratio 

3.3-20 

p 

Density of material 

3.3-4 

a 

Nodal stress vector 

3.3-4 


Yield stress 

3.3-4 


Stress tensor 

3.3-4 

Q 

Denotes a deformable body 

3.3-4 

3J2 1 

Denotes surface with prescribed displacements 

3.3-4 

3S2 2 

Denotes surface with prescribed tractions 

3,3-4 

Superscript 

Description 

Page 

n 

Iteration number 

3,3-6 

PL 

Denotes plastic quantity 

3.3-29 

T 

Denotes transpose of a matrix 

3,3-5 

1 

Denotes element coordinates 

3,3-9 

A 

Denotes a prescribed quantity 

3,3-4 

* 

Denotes a virtual quantity 

3,3-4 
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3.4 BOUNDARY ELEMENT METHOD 


3.4.1 Overview 

It is the goal of the Advanced Formulation Development portion of this program 
to develop a computational technique for the solution of linear, nonlinear and 
transient problems in gas turbine engine hot section components. This tech- 
nique is to be distinct from, and complementary to, the finite element method. 
The existence of such a computational tool will enhance the ability to cali- 
brate the other codes developed under this contract. In addition, it is to be 
expected that different techniques will prove optimal, in terms of efficiency 
or accuracy, for particular types of component analyses. Since almost all gen- 
eral purpose structural analysis computer programs presently available employ 
the displacement finite element method, the new program developed as part of 
the Advanced Formulation Development effort can be expected to extend the 
ability to perform realistic analyses of hot section components. 

Pratt & Whitney and its subcontractor, the State University of New York at 
Buffalo, have chosen, with the agreement of the NASA-Lewis program manager, to 
develop the boundary element method to fulfill the requirements of Task IC, 
and expect to continue its development throughout the remainder of the Inelas- 
tic Analysis Methods program. Very significant progress has been achieved dur- 
ing the first year of the contract effort. The analytical basis for a general 
purpose structural analysis program employing the boundary element method has 
been developed. Elastic, inelastic and both periodic and aperiodic dynamic ef- 
fects have been considered. Numerical methods for the implementation of these 
analyses have been selected or, in many cases, newly developed. A computer 
program incorporating these techniques with sufficient generality to allow the 
analysis of hot section components has been designed, developed and tested. 

The first year development effort is discussed in more detail below. Section 

3.4.2 reviews the results of a boundary element method literature survey con- 
ducted early in the program. Section 3.4.3 summarizes the analytical and nu- 
merical basis of the boundary element method for elastic, inelastic and dynam- 
ic problems in three dimensions. The overall structure and development of the 
boundary element method computer program is described in Section 3.4.4. Vali- 
dation/verification of the resulting code is reviewed in Section 3.4.5. 

3.4.2 Literature Survey 

3. 4. 2.1 Introduction 

This section contains a brief review of available literature on the boundary 
element method, applied to problems of linear, nonlinear and dynamic stress 
analysis. Attention is primarily focused here on the boundary element litera- 
ture relevant to the present contract. 
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The earliest work based on the direct boundary element formulation was pre- 
sented by Jaswon and Ponter (Reference 1) for the Saint-Venant torsion problem 
of elastic bars, and by Rizzo (Reference 2) for elastostatic problems. In re- 
cent years, Cruse, Rizzo, Lachat, Watson, Banerjee, Butterfield, Brebbia, 
Wilson and others have all contributed to development of the direct boundary 
element method. At present, higher-order elements, as in the finite element 
method, have been employed and some standardization of solution procedures has 
been achieved for linear problems in solid mechanics. A number of textbooks 
(References 3 through 7) and advanced level monographs (References 8 through 
10) provided not only a full description of these developments but also of de- 
velopments in other fields. Reports by Cruse and Wilson (References 11, 13 and 
14) on fracture mechanics, and Mendel son on plasticity (Reference 12) provide 
very readable accounts of developments in these areas. 

In the discussion below, the development of the boundary element method is 
outlined for several different fields of continuum mechanics. 

3. 4. 2. 2 Linear Stress Analysis 

The initial application of boundary element methods in elastostatics was to 
torsion problems of elastic bars. The pioneering work in this area was done by 
Jaswon and Ponter (Reference 1). Subsequent to their work, computational effi- 
ciency and accuracy have both been much improved. 

Extension of both direct and indirect integral equation methods to the static 
stress analysis of three-dimensional bulky solids followed shortly after the 
successful application of these methods to two-dimensional problems. Refer- 
ences 4 and 11 through 23 describe the work of a number of researchers. In 
particular, Lachat and Watson (Reference, 16) and Watson (Reference 17) de- 
scribe algorithms developed using the direct formulation applied to the stress 
analysis of thick cylindrical pressure vessels with holes. Banerjee (Refer- 
ences 18 and 19) and Banerjee and Butterfield (Reference 20) describe numeri- 
cal solutions to problems involving pile foundations as well as buried foun- 
dations of arbitrary shape. 

Cruse, Wilson and others (References 23 through 30) show excellent applica- 
tions of the boundary element method to two- and three-dimensional problems of 
linear elastic solids. Alarcon, Brebbia and Dominquez (Reference 31) showed an 
equivalence between the direct boundary element method and the method of re- 
siduals, and solved some simple problems. A more general discretization was 
presented by Lachat and Watson (References 16 and 17) by introducing higher- 
order boundary elements, whereby the efficiency and accuracy of numerical cal- 
culation were much improved. Anisotropic elastic problems were investigated by 
Wilson and Cruse (Reference 32). Elastic stress analysis of axi symmetric 
bodies was carried out by Cruse, Snow and Wilson (Reference 27), and also by 
Mayer, Drexler and Kuhn (Reference 33), under centrifugal body forces as well 
as steady state temperature gradients. 
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If the temperature distribution is known, the uncoupled thermoelastic problem 
can be treated as an elastostatic problem loaded by a body force corresponding 
to the temperature distribution. Rizzo and Shippy (Reference 34) used this 
approach to solve three-dimensional thermoelastic problems for nuclear pres- 
sure vessel components and for turbine airfoils. Chaudouet and Loubignac (Ref- 
erence 35) performed a three-dimensional analysis of the thermoelastic behav- 
ior of rollers. 

Hansen (Reference 36) and Stern (Reference 37) proposed direct boundary ele- 
ment methods for plate problems formulated in terms of variables having clear 
physical meanings. Independently, Bezine (References 38 through 40) also pro- 
posed similar methods of solution. Wu and Altiero (Reference 41) have recently 
proposed a new solution procedure for the elastic bending of anisotropic 
plates and solved some example problems. 

Stress intensity factors in elastic fracture mechanics have been calculated 
with considerable success using boundary element methods (References 23 
through 30, and 42 through 47). Results have been obtained for various crack- 
ing modes in both two and three dimensions. Recently, crack tip boundary ele- 
ments have been introduced which take into account the singular behavior of 
displacements and tractions at crack tips (References 14 through 47). A stress 
intensity factor was also computed by Snyder and Cruse (Reference 28) for ani- 
sotropic elastic plates. 

It is clear that solutions to problems in elastostatics based on boundary in- 
tegral equations have reached such a state of development that it is now pos- 
sible to undertake elastic and thermoelastic stress analyses of two- and 
three-dimensional bodies in a routine manner. For instance. Reference 4 dis- 
cusses the analysis, using an isoparametric boundary element representation, 
of a disc rim-slot which retains blades in a gas turbine rotor assembly. 

3. 4. 2. 3 Dynamic Stress Analysis 

As far as the solution of the general transient two-dimensional linear elasto- 
dynamic (or viscoelastodynamic) problem is concerned, there are currently 
three basic approaches available: 

1. Determination of the steady-state solution by boundary element method 
approaches followed by reconstitution of the transient response using 
Fourier synthesis, as done by Banaugh and Goldsmith (Reference 48), 
Kobayashi and Niwa, et al (References 49 and 50). 

2. Solution of the problem in the Laplace transform domain by boundary 
element method approaches followed by a numerical inverse transforma- 
tion to obtain the response in the time domain, as done by Doyle 
(Reference 51), Cruse and Rizzo (Reference 52), Cruse (Reference 53), 
and Manolis and Beskos (Reference 54). 
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3. Time domain formulation and solution of the problem by boundary ele- 
ment method approaches in conjunction with step-by-step time integra- 
tion schemes. This has been done by Cole, et al (Reference 56), for 
the anti -plane strain case and by Niwa, et al (Reference 55) for the 
two-dimensional case. 

A comparison of the above three approaches on the basis of their accuracy and 
efficiency was done by Manolis (Reference 57). It should be noted that in the 
above papers some simple two-dimensional elastodynamic problems were solved 
such as: 1) the case of an unlined or lined circular cylindrical cavity under 
the passage of longitudinal or transverse waves (References 49 and 50), 2) the 
cases of square or horseshoe shaped cylindrical cavities under longitudinal 
waves, and 3) the case of wave propagation in half-planes, etc. 

In conclusion, it may be mentioned that although the static version of the 
boundary element method is rather fully developed, this is not the case in 
elastodynamics. Although in the earlier attempts the dynamic verson of the 
boundary element method showed considerable promise, the development has not 
been carried out to the stage where the method could be considered as an 
established problem solving tool. 

3. 4, 2. 4 Nonlinear Stress Analysis 

Since material nonlinearities may be introduced into an elastic analysis meth- 
od as a set of initial strains, stresses or body forces in the same way that 
thermal loading or centrifugal effects can be introduced, it is hardly sur- 
prising that such a successful analysis tool for elastostatics can also be de- 
veloped for plasticity and creep analysis. However, for the boundary element 
method, unlike the finite element method, it is possible to provide a brief 
and yet quite comprehensive review of work published to date on nonlinear pro- 
blems. 

The earliest boundary element formulation for elastoplasticity was due to 
Swedlow and Cruse (Reference 58), Rzasnicki, et al (Reference 59) proposed an 
initial strain approach for the three-dimensional analysis of work hardening 
materials, but only planar problems were solved. 

Mendel son, Rzasnicki and Albers (Reference 59) also describe formulations for 
elastoplastic analysis of torsion and planar problems. Both a bi harmonic for- 
mulation based on the Airy stress function and the displacement (direct) for- 
mulation of Swedlow and Cruse were utilized in numerical solutions (generally 
for V-notched beams). The elastoplastic solution was of an initial strain type 
using the method of successive elastic solutions. Boundary and interior dis- 
cretizations were of the simplest type - straight boundary segments with the 
midpoint representing the boundary value on the element and an interior grid 
of essentially rectangles with strains evaluated at midpoints using the (anal- 
ytically evaluated) integral equations for strain. Recently, Telles (Refer- 
ences 60 and 61) has adopted a similar but slightly more sophisticated ap- 
proach for a variety of problems involving strain hardening and perfect plas- 
ticity. In his solutions, a linear variation of boundary displacement and 
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tractions is adopted, with interior triangular cells in which strains are com- 
puted at the vertices. Thus, a linear variation of strain across an interior 
cell is computed, with an associated improvement in the accuracy of modeling 
the plastic strain distribution. 

Mukherjee and coworkers (References 7, and 64 through 67) have utilized a rate 
formulation of the inelastic equations and a scheme for integrating forward in 
time. Material behavior was described either by the equations of power law 
creep or by a state variable consitutive model which draws no distinction be- 
tween classical time-independent plasticity and time-dependent creep. The an- 
elastic strain rates at t=0 are obtained from the consitutive equations, and a 
knowledge of these enables boundary displacement and traction rates and inter- 
ior displacement and stress rates to be obtained. The values of these quanti- 
ties at a new time t+At are then obtained by integrating forward in time using 
an Euler-type strategy with automatic time step control. In Reference 67, a 
comparison between boundary element and finite element solution times is made, 
with the boundary element program giving a faster and more accurate solution. 

Classical time-independent plasticity solutions have been successfully ob- 
tained by Banerjee, et al (References 68 through 72) using an incremental ini- 
tial stress procedure similar to that used for finite element analysis. Plane, 
axisymmetric and three-dimensional problems have been solved. The initial 
stress algorithm is described in some detail in the boundary element method 
book by Banerjee and Butterfield (Reference 8). This work differs in one par- 
ticular respect from that of other workers in that strains and stresses are 
calculated from displacements at the interior cell nodes rather than using the 
integral equation for strain. This introduces some loss in accuracy but is 
considerably more efficient computationally. In later work, boundary geometry 
and unknowns are represented by quadratic elements (References 71 and 75). 
Cathie and Banerjee (Reference 72) described a generalized formulation to take 
account of time-dependent inelastic deformation, conventional plasticity, vis- 
coplasticity and creep with a single inelastic algorithm. Other conventional 
elastoplasticity solutions have also been reported (References 73 and 74). 

The boundary element method has recently been extended to the analysis of 
large deformation problems (References 76 and 77). Although these publications 
essentially deal with very small test problems involving simple rectangular 
two-dimensional regions and do not develop the necessary surface or interior 
identities, feasibility has already been demonstrated. 
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3,4.3 Formulation Development 

3. 4. 3.1 Summary 

Various boundary element formulations for linear, nonlinear and dynamic pro- 
blems are outlined in this section. The major governing equations are defined 
and the relevant integral formulations are presented. 

An advanced formulation of the boundary element method has been developed for 
three-dimensional problems of elastoplasticity and viscoplasticity and for dy- 
namics. In this formulation both the surface geometry and unknowns have been 
represented by isoparametric boundary elements. An efficient mapping scheme 
has been devised such that the kernel-shape function-jacobian products can be 
efficiently and accurately evaluated using the standard Gaussian integration 
formulae. For nonlinear analysis, a fast and accurate solution algorithm has 
been developed where the previous history of the development of initial 
stresses and plastic strains is utilized. The dynamic analysis has been de- 
veloped in such a manner that the steady state (periodic) dynamic as well as 
transient forced vibration analysis can be carried out using the same imple- 
mentati on , 

3. 4. 3. 2 Linear and Nonlinear Stress Analysis 
Basic Governing Equations 


The governing differential equation for a solid in which the inelastic strain 

rate 1^. (i.e., plastic + thermal), and body force excitation rate f is 
present, can be expressed as: 
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where ii is the displacement rate 
e?. = plastic strain 

* V 


temperature 


f = body force per unit volume and 


dots on quantities represent their time derivatives. 
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The corresponding boundary integral formulation for the problem can be easily 
shown to be (see Banerjee and Butterfield, Reference 4): 
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For a linear stress analysis problem, the quantities within the volume 

integral sign, i.e., e^. (= {.. EaT) and f are both known. If the body forces 

fi“ are conservative, then the volume integral can, if desired, be converted 
into a surface integral. For nonlinear analysis, is also present, and must 
be determined at every step of the solution process. The plastic strain is 
usually a function of the state variables such as the stresses cr-jj, plastic 
strains e?. and temperature T, 


Although equation (3.4-2) is written for a surface point Iq* '''t can also be 
used for an interior point | if the term C^-j arising from the treatment of 
the improper integral involving the function F^j is assumed to be zero, i.e.. 
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(3.4-3) 
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UjU) = ^Ct.(x)G^.j(x,|) - F^.j(x J)u^.(x)] dS 


+ /T.|^j.(x,ne5k(x)dV + ^G.j.(xJ)f.(x)dV 

By calculating the strains corresponding to equation (3.4-3) and substituting 
In the stress-strain relations: 


”jk = «jk - (i§|i (3.4-4) 


2uv 


we can obtain the following expressions for the stresses at an Interior point: 
^jk(^) " sft^-(x)D..jl^(x,^) - S^.j|^(x J)u^.(x)] dS 
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where D Is a small sphere around the source point i and the free term on the 
right hand side Is the result of an exact Integration over the volume D. 
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If the strains are purely devlatoric, then = 0 and equation (3.4-5) can be 
simplified, resulting In a small saving In computation. 


Equation (3.4-5) Is not only Invalid on the surface but also Is difficult to 
evaluate numerically at points close to It, For points on the surface, the 
stresses can be calculated by constructing a local Cartesian coordinate system 
with the axes 1 and 2 directed along the tangential directions and the axis 3 
In the direction of the outward normal. The stresses referred to these local 
axes (Indicated by overbars) are then given by: 


3.4-8 



'11 


1-v "^3 , 2 ^^11 ^22^ 1+v ^11 

1-v 


°12 = °Z1 = ITT+Vr ^12 


V X 


Ev 


'22 = ^3 ^ 7 "“^ ^"11 ^ " 22 ^ ^ "22 

1-v 


''33 = "^3 


032 - 023 = ^2 


''31 = ''13 = h 


(3.4-6) 


where E is the Young's modulus and I-jj defines the components of the elastic 
strains (i.e., total -inelastic) in the local axes system. This method of eval- 
uating the stresses on the surface was originally devised by Cruse (Reference 
25). 

The development presented here has assumed that all of the integral equations 
discussed are written over the entire structure to be analyzed. This restric- 
tion is not essential. In fact, in practical analysis it is almost always de- 
sirable to subdivide a complex structure into two or more subregions (referred 
to in this program as generic Modeling Regions, or GMRs). The integral equa- 
tions are written separately on each GMR, and the overall structure is then 
assembled by enforcing appropriate compatibility and equilibrium conditions on 
the boundaries between GMRs. This tactic can improve the numerical stability 
of the resulting equation system and also results in very considerable re- 
ductions in computer time for both elastic and inelastic problems. 

Solution Algorithm 

If we discretize the boundary using n isoparametric quadratic surface elements 
and the likely elastoplastic region using m isoparametric quadratic hexahedral 
cells, then the complete state of the structure can be defined in terms of the 
nodal values of displacement and traction. The surface and volume integrals of 
equation (3.4-2) can then be expressed as a sum of integrals over surface ele- 
ments and volume cells. Each of these integrals can be evaluated numerically. 
We can thus derive an algebraic representation of equation (3.4-2) as: 
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A X = B y + C 


(3.4-7) 


where x and y are unknown and known boundary tractions and displacements, 
respectively, and: 

C = G f + T 

couples the effect of the volume cells to the boundary solution. 

It should be emphasized that for a linear analysis, C is completely prescribed 
and equation (3,4-7) can be solved directly. ~ 

Similarly, the interior equations for displacements and stresses may be cast 
in the form: 


u^ = X 


y + 


(3.4-8) 
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where the superscript i is used to designate interior field point. 

The solutions for a nonlinear analysis can be developed as an accelerated 
elastoplastic algorithm or a viscoplastic algorithm depending on the nature of 
the constitutive model involved. For conventional elastoplastic constitutive 


relations, an initial stress (defined as 5°^ = where 


is the elastic 


constitutive matrix) algorithm would be as follows: 

1. Obtain an elastic solution and scale to first yield. Hence, determine 


.1 

O 


* * i 1 

2. Evaluate x, u , a for a small load increment with an initially 

estimated value of a° obtained from the extrapolation of the 
previously generated history of initial stresses (if no prior plastic 
history exists, assume a® = 0) 
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3. Accumulate all incremental quantities of stress, displacement 

4. Calculate with the new stress history the current constitutive matrix 
and determine the initial stresses 


11 

5. Calculate new x, u , a with y = 0 (i.e., no boundary loading change) 


6. Return to step 3 if x is greater than an acceptable norm 

7. Return to step 2, the next load increment. 

During the early part of the loading, the algorithm will require a few itera- 
tions for each load increment {since no history of generation of initial 
stresses exists). From the second or third increment onwards the solution con- 
verges after one or two iterations. 

If a viscoplastic model, including thermally induced plastic creep strain be- 
havior, is used the solution is best obtained using a current rather than an 
incremental form of equations (3.4-7) through (3.4-9). These equations, if 
written in non-rate form, are valid expressions of the elastic problem given 
the boundary conditions and an initial strain field. Providing the developing 
initial strain field in the plasticity and creep problem is integrated in a 
sufficiently accurate manner, the solution to the nonlinear problem at any 
time may be obtained using current values. Advantages are that the algorithm 
is simplified, errors introduced through the accumulation of incremental quan- 
tities are minimized and the equilibrium equations are satisfied at each stage 
of the calculation. 

The essentials of the viscoplastic solution algorithm are as follows: 

1. Evaluate x = A“^ B y + C 

o'' = X + s’ y + m’ 

for the prescribed boundary loading with the vector y and initial 
viscoplastic strain = 0 


2. Scale elastic solution to yield point for plasticity problem and 

apply a small load increment xy where x is fraction of the total load 

to be applied in an incremental analysis 
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3. Evaluate using the constitutive equation at the time t 

4. Select At and compute from a suitable integration of at time 
t + At 

5. Resolve equations in step 1 

6. Return to step 3 until is negligible or t > t„,„, the time up to 

luaX 

which the solution is required. 

5. For plasticity applications return to step 2 for the next load 
increment. 

The algorithm outlined above would be successful if the chosen 

where ^^cy^tical obtained from the careful study of the governing stiff 

differential equation; 


vp 

= -^ = F(a,eP) (3.4-10) 

A good working rule often adopted in practice is that the creep strain rate 
must not exceed a certain fraction of the total elastic strains at a critical 
section of the structure, i.e.. 





a U o 


where C® is 


the elastic compliance matrix 


a is the stress state 


a is a value between 0.1 to 0.5. 
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3. 4. 3. 3 Dynamic Stress Analysis 
Governing Equations 


The governing differential equations for the response of a three-dimensional 
solid in the absence of body forces can be expressed as: 


(cl - c\) u.^.j + c\ Uj^.. = u. 


(3.4-12) 


where Uj(xi‘,t) are the components of the displacement vector, indices i 
and j correspond to Cartesian coordinates x-f, commas indicate differentiation 
with respect to these coordinates, dots indicate differentiation with respect 
to time t and summation over repeated indices is understood; ci and C 2 are 
the propagation velocities of dilatational (P-wave) and distorsional (S-wave) 
waves, respectively, given in terms of the Lame constants x and u and the mass 
density p of the material by 


Cj = (x + 2ii) / p C 2 = p/p 


(3.4-13) 


The constitutive equations are of the form 


a • • s 0 


(cl - 2c 2^ ‘^r,rS‘j ^2 (“i,j “j,i) (3.4-14) 

where 0 'fj(xT,t) are components of the stress tensor and is the 

Kronecker delta. 


The Laplace transform with respect to time of a function f(xi,t) is defined 
as: 


L(f) = f(x.,s) = / f(x.,t)e"^^dt (3.4-15) 

where s is the transform parameter. 


Application of 
(3.4-14) under 

the 

zero 

Laplace transform to equations (3.4-12), 
initial conditions yields 

(3.4-13) and 



- 4 

- 2 - 2- 

i u. ..+CoU. ..-s u.=0 

7 1.10 2 j,n j 

(3.4-16) 


^ij 

-.[{4 

- ^'^0 Vr«iJ * 

(3.4-17) 





(3.4-18) 


3 . 4-13 



If a homogeneous, isotropic and linear elastic body is under the influence of 
harmonic forcing functions of time prescribed through boundary conditions of 
the form 

t.{K^,t) = e^'“^ (3.4-19) 


where u is the circular frequency and t^ the amplitude of the forces, then 
the displacement vector will be of the form 


Uj(X|^,t) = Qj(x,^,u) e^*“* 


(3.4-20) 


where u-,* is its amplitude. Substitution of equation (3.4-20) in (3,4-12) and 
(3,4-14) and cancellation of the common factor e''“t yields 
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where o’fj denotes the stress amplitude. A comparison of equations (3.4-16) - 
(3,4-18) with equations (3.4-21) - (3,4-23) indicates that any method for 
solving the Laplace transformed general transient problem can also be used to 
solve the steady-state problem as a special case, if the complex Laplace 
transform parameter s is replaced by -iu in the formulation. The solution u-f 
and a-ij will be, of course, in terms of the frequency u so that no inversion 
is required. The approach implies that the formulation corresponding to algo- 
rithms of Laplace inversion working with complex data will be used. 


The corresponding boundary element formulation in the transformed space 
(s,Xi*) follows from the above differential equation. Thus for a boundary 
point ^o» we have 
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dS (3.4-24) 


where s is either the transform parameter or -iu, and u is the circular fre- 
quency 
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Note that although the functions G and F become identical to their static 
counterparts G and F as u tends to zero It is important to evaluate this limit 
carefully because of the presence of u in denominator. As u increases these 
functions remain well behaved. 

Equation (3.4-24) can also be written for an interior point ^ by simply assum- 
ing Cjj = 0 and ^ = ^o- This gives the required identity for calculating 
the interior displacements. The interior stresses can be obtained from: 
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The corresponding equations for the stresses on the surface are identical with 
those for the static case with dynamic elastic constants given by equation 
(3.4-13). 


Numerical Implementation 


Since the basic governing equations for a dynamic analysis in the transformed 
space (either in s or u space) are virtually identical to the corresponding 
equations for the static analysis, the numerical implementation developed for 
the static case can be used to extract solutions for the dynamic problem for 
one value of the transform parameter s or frequency parameter u. It should 
also be noted that any internal viscous dissipation of energy (damping) can 
easily be accounted for by replacing the elastic moduli \ and 4 by their 
complex counterparts x* and u*: 
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U* = u(l + 2 ib) 

X* = x{l + 2 ib) 


(3.4-26) 


leaving Poisson's ratio unaltered. 

The parameter 2s, by analogy with single degree of^ freedom systems may be 
termed the damping ratio and is given by: 


B = (Oti/2u 


(3.4-27) 


where n is the coefficient of viscosity in the Kelvin-Voigt model. 

For a general transient case the solution must be transformed back to real 
time by a suitable inversion process. Let ^(t) be a real function of t, with 
f(t) = 0 for t < 0. The Laplace transform f(s) of f(t) is defined by equation 
(3.4-15) and its inversion formula is given by: 


X+ioo 




St 


ds 


V-i 


X-loo 


(3.4-28) 


where x(>0) is arbitrary, but greater than the real part of all the singular- 
ities of f(s), and s is a complex number with Re(s)B > 0. 

In the present applications f(s) is too complicated to be inverted analytic- 
ally and is available in numerical form. In those cases a numerical inversion 
of the Laplace transform is imperative. 


The method of inversion used is that of Durbin (Reference 78) and is actually 
an efficient improvement of the method of Dubner and Abate (Reference 79) 
which is based on the finite Fourier cosine transform. Durbin (Reference 78) 
combined both finite Fourier cosine and sine transforms to obtain the inver- 
sion formula: 


f(tj) 


where 




M-1 

4 Re[f(6)] + Re ( S (A(n) + iB(n))W 
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(3.4-29) 


S„ = 6 + in ~ 1 = V-1 W = e 


i2ir/N 
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B( 


n) = 2 Re{f(B + i(n + £N) ^)} 
1=0 

n) = S Im|f^B + i(n + £N) 


(3.4-30) 


n = 0, 1, 2, 3,...,N £= 0, 1, 2, 3,...,L 
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and f(t) is computed for N equidistant points tj = jAt = jT/N,j =0, 1, 2 , 
It is suggested that for LxN ranging from 50 to 5000 one should se- 
lect 6l = 5 to 10 for good results, where T is the total time interval of in- 
terest. The computations involved in equation (3.4-29) are performed by em- 
ploying the Fast Fourier algorithm of Cooley and Tukey (Reference 80). The al- 
gorithm employed here was used with considerable success by Manolis and Beskos 
(Reference 54). In the case of the Fourier transform (s = -iu), the above al- 
gorithm is equivalent to a Fourier synthesis. In this work we have used L=l, 
N=50 and eT=6, although as mentioned above these values can be adjusted to im- 
prove accuracy. 

The only topic remaining for discussion is the direct transform, i.e., to find 
the Laplace or Fourier transform of a given forcing function in time. If the 
forcing function is piecewise linear in time, which is an excellent approxima- 
tion for functions that are densely sampled, then the following exact formula 
can be used: 


f(s) 


N-1 

Z 

i=l 





+ S At 
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where F-,- is the value of f at time t-f and aF = F^+i - F-f. 

It should be noted that the Laplace transform solution is essentially a super- 
position of a series of steady state solutions and is therefore applicable 
only to linear elasto-dynamic problems. For nonlinear problems the solution 
must be obtained in the real time domain. The boundary integral formulation 
for the time embedded elasto-dynamic analysis has already been developed. This 
formulation can be extended to deal with nonlinear dynamic problems. 

The present steady state dynamic algorithm can also be used to determine the 
various vibration modes of a structural component. Unfortunately, the algo- 
rithm is computational ly inefficient. A new real variable version of the 
eigenvalue analysis is currently being developed which is likely to prove 
useful. The present steady state dynamic analysis is ideally suited for the 
stress analysis of components at a given value of frequency parameter u. 

For buckling analysis, the necessary boundary integral formulation has not yet 
been developed. Although it is possible to construct an integral representa- 
tion, it is doubtful if it can be made sufficiently computationally efficient 
to be a viable alternative. 
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3.4.4 Computer Program Development 
3. 4. 4.1 Introduction 


The mathematical goal of the computer program development is the accurate and 
efficient implementation of the analysis described in Section 3.4.3. By it- 
self, this is a significant task. While many of the capabilities required have 
been demonstrated in special purpose computer programs or discussed in the 
open literature, no program (even for two-dimensional analysis) includes all 
of the required capabilities. The construction of a three-dimensional code in- 
cluding elastic, nonlinear and dynamic capabilities was a most challenging 
task. 

Of equal importance is the degree of generality required in the definition of 
component geometry, loading and material properties. This is necessary if the 
program is to be applicable to real problems in the aerospace industry. Since 
no general purpose boundary element program existed as a starting point, or 
even a model, it was necessary to develop, during a single year, a new general 
purpose system for structural analysis. 

During the past year, the computer program BEST (J|oundary Hement ^tress Tech- 
nology) has been developed. It is capable of carrying out elastic, inelastic 
and dynamic analyses for real component geometries. The capabilities of the 
program have been calibrated by comparison with analytical solutions, other 
numerical stress analyses and data. In addition, the code was written with 
sufficient generality to allow the incorporation of new technology, developed 
during the remainder of this contract, without major recoding of the already 
completed portions of the program. 

The development of the computer program BEST is discussed in the following 
sections. 


3, 4. 4. 2 Global Program Structure 

BEST is designed to be a fully general structural analysis system employing 
the boundary element method. The program is written using standard IBM FORTRAN 
IV. Development has been carried out at Pratt & Whitney using an IBM 3081 in 
both batch and interactive modes and at SUNY-B on an HP9000 minicomputer sys- 
tem. In both cases the required code and workspace fit in core without re- 
quirement for overlays. The nature of the method is such that, for any real- 
istic problem, not all required data can reside simultaneously in core. For 
this reason extensive use is made of both sequential and direct access scratch 
files. 

The overall structure of BEST is shown in Figure 3.4-1. The program first ex- 
ecutes an input segment which is common to all types of analysis. After the 
input has been processed, there are three major branches, corresponding to 
elastic, inelastic and dynamic analysis. 
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INPUT PROCESSING 


STATIC ^X^TATIC0R\ DYNAMIC (COMPLEX VARIABLE CODE) 

\DYNAMIC?^^ I 


SURFACE INTEGRATION 


SURFACE INTEGRATION 


VOLUME INTEGRATION 
(IF REQUIRED) 


TRANSFORM BOUNDARY 
CONDITIONS 
(IF REQUIRED) 


SYSTEM EQUATION ASSEMBLY 


SYSTEM EQUATION ASSEMBLY 


SYSTEM MATRIX DECOMPOSITION 


SYSTEM MATRIX DECOMPOSITION 




Figure 3.4-1 BEST - Overall Program Structure 
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If only an elastic analysis is required, the necessary surface integrals are 
calculated and assembled into the set of system equations using specified 
boundary conditions. The system matrix is then decomposed and saved on disk, 
followed by the calculation of the solution vector. The full displacement and 
traction solution on each boundary element is then reconstructed from the so- 
lution vector. In a time dependent problem the process of constructing the 
load vector for the system equations is repeated at each time step, but the 
integration, formation and decomposition of the system matrix are done only 
once. 

If a dynamic analysis is to be carried out, a branch is taken immediately 
after the input processing. The main differences between the elastic and dy- 
namic portions of the code are: 

1. in the calculation of Laplace transforms of the boundary conditions 
(in the aperiodic case), 

2. in the use of complex versions of the quasi -static elastic routines 
connected with surface integration, system matrix assembly and system 
matrix decomposition, and 

3. in the use of the transform of the (complex) dynamic point load solu- 
tion in place of the Kelvin solution. 

In the case of inelastic analysis the elastic branch is followed through the 
calculation of the surface integrals required for both system matrix creation 
and the evaluation of elastic stresses at the points used to describe the vol- 
ume distribution of inelastic strain. The system matrix is decomposed once, 
for the initial elastic solution, and the decomposed form is saved for re- 
peated use during the plasticity calculations. The main difference between the 
elastic and inelastic analyses is the iterative process (carried out at each 
load or time step) in which the inelastic strain distribution is evaluated and 
employed in volume integrals which modify the system equation load vector. 
This process is nested inside the time stepping procedure of the elastic code. 

It should be noted that, in the present version of BEST, the inelastic volume 
integration routines are used to deal with both thermal stresses (due to 
either transient or steady state temperature fields) and inhomogeneous elastic 
material properties. Effort during the second year of the program will be de- 
voted to developing a more efficient treatment of these effects in purely 
elastic problems. 

Various aspects of the computer program are discussed below. 
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3. 4. 4. 3 Program Input 

The input for BEST is free field. Meaningful keywords are used to identify 
data types and to name particular data sets. The overall input structure is 
shovi/n in Figure 3.4-2. The input is divided into four major types: 

1 . Case Control Cards 

The case control cards define global characteristics of the problem. In ad- 
dition to the problem title, these input items define the ways in which the 
problem departs from a static, homogeneous elastic analysis. The reading or 
writing of restart data is also defined at this point. In all cases the ab- 
sence of a case control card will cause the program to default to the simplest 
case, a standard elastic analysis. 

2. Material Property Definition 

The material property input allows the definition of elastic and, if required, 
inelastic material properties for a variety of materials. For isotropic, homo- 
geneous materials both Young's modulus and the coefficient of thermal expan- 
sion can be prescribed in tabular form for a user-defined set of temperatures. 
Temperature independent values of density and Poisson's ratio are defined. 

If inelastic material properties are required, the particular constitutive 
model to be used must first be chosen. The constants required to initialize 
the model for a given material are then input, together with a material iden- 
tifier. 

For both elastic and inelastic properties, provision is made to link BEST to 
existing user material properties libraries. In this case the user is respon- 
sible for providing code to access the local libraries, 

3. Geometry Input 

Geometry input is defined one GMR (generic modeling region, or subregion) at a 
time. To initiate the input, a tag is provided to identify the GMR, a material 
name and reference temperature are defined to allow initialization of material 
properties and various flags are input (if required) to identify behavior 
other than static elastic within the GMR. 

The next block of geometry input consists of the Cartesian coordinates of the 
user input points for surface and volume geometry definition, together with 
identifiers (normally positive integers) for these geometric nodes. 
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Following the definition of an initial set of nodal points, the surface con- 
nectivity of the GMR is defined through the input of one or more named sur- 
faces. Each surface is made up of a number of elements, with each element de- 
fined in terms of several geometric nodes. All surface element types presently 
employed represent surface geometry using the quadratic isoparametric shape 
functions (Figure 3.4-3). Three sided elements, defined using six rather than 
eight geometric nodes, are used for mesh 'transition purposes. The terms quad- 
rilateral and triangle are normally used to refer to the eight and six noded 
elements, although the real geometry represented is, in general, a nonplanar 
surface patch. 



NO 


INTERFACE 

INPUT 


NO 


BODY FORCE 
INPUT 



1 



END 





Figure 3.4-2 BEST - Input Structure 
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Figure 3.4-3 Quadratic Isoparametric Shape Functions 


Over each element the variation of displacement and traction can be defined 
using either the linear or quadratic isoparametric shape functions. Linear and 
quadratic elements can share a common side, which is then constrained to have 
linear displacement and traction variation. 

In addition to the element types mentioned above, elements which extend to in- 
finity are provided. These elements are designed to allow modeling of struc- 
tures connected to ground, and automatically incorporate appropriate decay 
conditions. The characteristics of the various element types are summarized 
bel ow. 


Element Type 

Geometry Nodes 

Displacement/Traction Nodes 

Linear Quadrilateral 

8 

4 

Linear Triangle 

6 

3 

Quadratic Quadrilateral 

8 

8 

Quadratic Triangle 

6 

6 

Linear Infinite 

8 

2 

Quadratic Infinite 

8 

3 


The geometry input includes the option to duplicate surfaces. Once a given 
surface has been completely defined, a copy of it can be automatically cre- 
ated, with arbitrary translation and rotation. In this case the program will 
internally define and number any new nodes required, and will construct all 
element definitions. The surface referenced can be from the current GMR, or 
from any GMR previously defined. 
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The final input for each GMR is the definition of volume cells (presently re- 
stricted to twenty node isoparametric elements) for that portion of the GMR in 
which volume effects (inelastic or thermal) are anticipated. The nodes of the 
(partial) volume discretization lying on the surface of the GMR need not be 
coincident with nodes of the surface discretization. 

4. Boundary Condition Input 

The final input section provides for the definition of boundary conditions, as 
functions of both position and time,. Data can be input for an entire surface, 
or for a subset (elements or nodes) of a surface. Input can be in global coor- 
dinates, or can define rollers or pressure in the local coordinate system. In- 
put simplifications are available for the frequently occurring cases of bound- 
ary data which is constant with respect to space and/or time variation. Each 
boundary condition set can be defined at a different set of times. 

Special types of boundary conditions which are available include interfaces 
(fixed or sliding contact) between two GMRs, cyclic surfaces (to allow effi- 
cient modeling of structures with periodic geometry and loading) and springs 
to ground. 

In addition to the boundary condition types defined above, two types of body 
force loading (thermal and centrifugal) can be defined. 

3. 4. 4. 4 Surface Integral Calculation 

Following the processing of the input data, the surface integrals occurring in 
equations (3.4-2) and (3.4-3) are evaluated numerically. This is the most time 
consuming portion of most elastic analyses, and contributes heavily to the 
cost of inelastic analysis. In BEST the results of these integrations are 
stored as they are calculated, rather than being assembled into the final 
equation system immediately. Although this is somewhat more costly in terms of 
storage and CPU (central processing unit) time, it has led to much greater 
clarity in the writing of the initial version of BEST. In addition, it pro- 
vides much greater flexibility in the implementation of various restart and 
boundary condition options. 

The basic outline of the surface integration process is shown in Figure 3.4-4. 
The calculations proceed first by GMR (generic modeling region), then by 
source point (the equation being constructed) and finally by surface element. 
The results for each source point/surface element pair are written to disk. 
All of the calculations are carried out and stored in the global (Cartesian) 
coordinate system. 

The calculations connected with the creation of the system equations are the 
more complex. In this case either singular or nonsingular integrals can be en- 
countered. The integrals are singular if the source point for the equations 
being con^structed lies on the element being integrated. Otherwise, the inte- 
grals are nonsingular, although numerical evaluation is still difficult if the 
source point and the element being integrated are close together. 
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Figure 3.4-4 BEST - Surface Integral Calculations 
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In both the singular and nonsingular cases Gaussian integration is used. The 
basic technique is developed in Reference 4, and was first applied in the 
three-dimensional boundary element method by Lachat and Watson (Reference 16). 
In the nonsingular case an approximate error estimate for the integral was de- 
veloped based on the work of Stroud and Secrest (Reference 81). This allows 
the determination of element subdivisions and orders of Gaussian integration 
which will retain a consistent level of error throughout the structure. Numer- 
ical tests have shown that the use of 3x3, 4x4 and 5x5 Gauss rules provide the 
best combination of accuracy and efficiency. In the present code the 4x4 rule 
is used for nonsingular integration, and error is controlled through element 
subdivision. The origin of the element subdivision is taken to be the closest 
point to the source point on the element being integrated. 

If the source point is very close to the element being integrated, the use of 
a uniform subdivision of the element can lead to excessive computing time. 
This frequently happens in the case of aerospace structures, due either to 
mesh transitions or to the analysis of thin walled structures. In order to im- 
prove efficiency, while retaining accuracy, a graded element subdivision was 
employed. Based on one-dimensional tests, it was found that the subelement 
divisions could be allowed to grow geometrically away from the origin of the 
element subdivision. Numerical tests on a complex three-dimensional problem 
have shown that a mesh expansion factor as high as 4.0 can be employed without 
significant degradation of accuracy. 

In the case of singular integration (source point on the element being inte- 
grated) the element is first divided into triangular subelements. The integra- 
tion over each subelement is carried out in a polar coordinate system with 
origin at the source point. This coordinate transformation produces nonsin- 
gular behavior in all except one of the required integrals. Normal Gauss rules 
can then be employed. The remaining integral (that of the traction kernel 
Tij times the isoparametric shape function which is 1.0 at the source point) 
is still singular, and cannot be numerically evaluated with reasonable effi- 
ciency and accuracy. Its calculation is carried out indirectly, using the fact 
that the stresses due to a rigid body translation are zero (Reference 16). It 
has been found that subdivision in the circumferential direction can be re- 
quired to preserve accuracy in the singular integration. A maximum included 
angle of 15 degrees is used. Subdivision in the radial direction has not been 
required. 

The modifications required in the surface integrals for the solution of dynam- 
ic problems are primarily the replacement of a number of routines with complex 
variable versions of the same code and the substitution of the transform of 
the dynamic kernel function. Modification of the calculation of the singular T 
term is also required, since the rigid body argument cannot be used directly. 
It has been found, however, that the singular behavior of this term is en- 
tirely due to the static kernel. The modifications required for the dynamic 
solution can be calculated entirely in terms of nonsingular numerical inte- 
grals. 
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The surface integrals required for calculation of displacement and stress at 
interior or surface points are of the same type as those involved in the gen- 
eration of the system equations, except that only nonsingular integrals are 
involved. If the source point involved is located on the surface of the body, 
then numerical integration is not required. Instead, the required quantities 
are calculated using the displacements and tractions on the element (or ele- 
ments) containing the source point, as discussed in Section 3. 4. 3. 2. 

3. 4. 4. 5 Volume Integral Calculations 

The calculation of volume integrals is required in both the construction of 
the system equations and the evaluation of interior stresses and displace- 
ments. In a purely elastic context this requirement is due primarily to the 
presence of body forces. In the case of inelastic analysis the inelastic 
strain distribution must also be integrated for each source point/volume cell 
pair. 

The three volume integrals involved in the analysis: G, T and M, all exhibit 
singular behavior as the load point and the source point coincide. Thfi order 
of the singularities for G, T and M are (1/r), (1/r^) and (l/r^), re- 
spectively. 

The overall strategy used in the volume integration is very similar to that 
described for the surface integration in the preceding section; i.e., based on 
the consistent level of error throughout the volume, Gaussian integration has 
been used to evaluate both the singular and nonsingular integrals. 

Nonsingular integrals involving T and M have been evaluated using 3x3x3, 4x4x4 
or 5x5x5 Gauss rules depending on the distance between the source point and 
field point. If the source point lies in the neighboring cell, the cell is 
subdivided into several sub-cells. The corresponding integrals involving the 
weakly singular function G have been evaluated using Gaussian integration 
rules of one order less than those mentioned above. 

Singular integrals involving the functions G and T were evaluated using a 
spherical polar mapping with the origin at the field point. The cell is sub- 
divided into several tetrahedra through the field point and the resulting 
sub-cell integrations are carried out by mapping each of these sub-cells into 
a unit cube. The kernel -shape function products are thus made to behave like a 
regular integrand. Gauss rules of 3x3x3 (for G) and 4x4x4 (for T) are used for 
evaluating these sub-cell integrals. 

Singular integrals involving the functions M times the shape function N re- 
quire some very special attention. For node at which the shape function is 
zero, the behavior of these integrands is similar to the singular integrals 
involving the function T as described previously; and accordingly they have 
been treated in a similar manner. Those involving the shape function attaining 
a value of 1,0 at the singular point (a situation which is equivalent to the 
case of traction kernel F times the shape function 1.0 at a singular node at 
the surface) cannot be evaluated accurately without a massive amount of com- 
putational effort. Fortunately, the free term obtained by the exact integra- 
tion over a spherical exclusion (see equation 3.4-5) accounts for 95% to 98% 
of the contribution from this integral; accordingly, this integral has been 
assumed to be zero. 
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3. 4. 4. 6 System Matrix Assembly 

The assembly of the system equations is a multi step process, as shown in Fig- 
ure 3.4-5. The process is the same for elastic, inelastic and dynamic pro- 
blems. The assembly process is carried out differently in the present program 
than in previous two- and three-dimensional codes. This is due to the fact 
that the present program is aimed largely at the solution of nonlinear pro- 
blems. Since the nonlinear effects are accounted for in the load terms of the 
equation system, and since these effects depend both on the present load state 
and on the past history of the system, it is not possible to calculate all of 
the system load vectors prior to the beginning of the solution process. As a 
result, the assembly process must be designed to facilitate the updating of 
the load vector, based on material response as well as time dependent loads. 



Figure 3.4-5 BEST - System Matrix Assembly 









The first step in the assembly process is the reduction of the rectangular ma- 
trix of T integrals to a square matrix. This matrix is the prototype of the 
system matrix. The columns of this matrix are transformed or replaced, as re- 
quired by the boundary conditions, as the assembly process proceeds. 

A key problem in the entire process is the proper definition of appropriate 
coordinate systems, on a nodal basis. This is a problem common to any direct 
boundary element method program which treats structures with nonsmooth sur- 
faces. It arises because the tractions at a point are not uniquely determined 
unless the normal direction to the surface varies continuously at the point in 
question. 

The original surface integral calculations are all done in global coordinates. 
If a displacement boundary condition is specified at a given node, in global 
coordinates, then no new coordinate system definition is required. It is only 
necessary to keep track of the subset of elements, containing the given node, 
on which the fixed displacement is to be reacted. However, if a displacement 
is specified in a nonglobal direction at a given node, then a new nodal coor- 
dinate system must be defined and, potentially, updated as further boundary 
conditions are processed. The associated nonzero reactions must then be ex- 
pressed in the new coordinate system. 

After all required local coordinate systems have been defined, any modified 
boundary conditions required by the presence of body forces are calculated. 
There are a variety of methods of accounting for body force distributions. The 
one employed in the multi -GMR version of BEST is a modification of the inte- 
gral equations resulting from the representation of the solution as the sum of 
a complementary solution and a particular integral. The boundary conditions 
must be modified to reflect displacements or tractions due to the particular 
integral. This results in changes to any explicitly specified boundary con- 
ditions and in the introduction of potentially nonzero boundary conditions on 
any surfaces previously defined implicitly as being traction free. 

Following this preparatory work, the final assembly of the system equations is 
carried out. It is performed in three major steps; 

1. Transformation of the columns of the matrices to appropriate local 
coordinate systems and incorporation of any boundary conditions in- 
volving springs. 

2. Incorporation of compatibility and equilibrium conditions on inter- 
faces between GMRs and on cyclic surfaces. On interfaces either a 
completely glued condition (full displacement compatibility) or a 
sliding condition (only normal displacement compatibility) is avail- 
able. 

3. Application of specified displacements and tractions. 
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Two particular features of the equation assembly deserve special comment. 
First, in multi-GMR problems the system matrix is not full. Rather, it can be 
thought of as consisting of an NxN array of submatrices, each of which is 
either fully populated or completely zero. Only the nonzero portions of the 
system equations are preserved during system matrix assembly. In order to im- 
prove the numerical conditioning of the system matrix for the solution pro- 
cess, the columns are reordered to number variables lying on the same inter- 
face, but belonging to two different GMRs, as close together as possible. The 
rows of the system matrix are placed in the same order as the columns. 

Second, rather than simply assembling an explicit load vector at each time 
point in the solution process, load vector coefficient matrices are assembled 
and stored. These allow the updating of the load vector at any required time 
point simply by interpolating the time dependent boundary conditions and per- 
forming a matrix multiplication. This capability is particularly important in 
inelastic analysis, since frequent updating of the load vector is required. A 
similar process is used in the calculation of interior and boundary stresses. 

3.4.4. 7 System Equation Solution 

A new solver was written for BEST. It operates at the submatrix level, using 
software from the UNPACK package (Reference 82) to carry out all operations 
on submatrices. The system matrix is stored, by submatrices, on a direct ac- 
cess file. The decomposition process is a Gaussian reduction to upper triangu- 
lar (submatrix) form. The row operations required during the decomposition are 
stored in the space originally occupied by the lower triangle of the system 
matrix. Pivoting of rows within diagonal submatrices is permitted. 

The calculation of the solution vector is carried out by a separate subrou- 
tine, using the decomposed form of the system matrix from the direct access 
file. The process of repeated solution, required for problems with time de- 
pendent and/or nonlinear behavior, is highly efficient. 

3. 4. 4. 8 Inelastic Solution Process 

The inelastic solution algorithm starts with an elastic analysis of the pro- 
blem for the first loading increment (complete with the specified boundary and 
body force loading). At the end of the elastic increment the state variables 
are calculated and the nonlinear consitutive equations are established. The 
difference between the correct stresses or strains and the elastic stresses or 
strains are estimated. These are used to compile a new boundary condition 
vector (the right-hand side of the system equation). The process is essen- 
tially repeated until the constitutive equations yield the stresses that are 
negligibly different from the calculated stresses using equation (3.4-5). A 
new loading increment is taken and the process is repeated for each subsequent 
increment (see Section 3. 4. 3. 2). 
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During such an incremental solution process, a considerable saving in the com- 
puting time has been achieved by introducing: 


a. a solver which allows for resolution of the same system equations for 
different right-hand sides so that the decomposition of the system 
matrix is done only once, 

b. a routine for the extrapolation of the plastic strains based on the 
previous loading history and the generation history of the plastic 
strains. 

Although the extrapolation of the plastic strains does sometimes affect the 
accuracy of the nonlinear solution, the savings could be very considerable 
since usually one iteration may lead to 95% of the desired results. 

The nonlinear solution algorithm and the associated routines are included 
within the program BEST in such a manner that their presence will not affect 
the efficiency of an elastic analysis. 

Various levels of constitutive relations have been programmed within BEST. 
These include a kinematic elastoplastic hardening model, and a viscoplastic 
kinematic hardening model of Walker admitting thermally dependent hardening. 

3. 4. 4. 9 Output 

The output from BEST is relatively straightforward. It consists of nine 
sections, as follows: 

1. Complete echo of the input data set. 

2. Summary of case control and material property input. 

3. Complete definition for each GMR, including all surface and volume 
nodes, surface elements and volume cells. 

4. Complete summary for each boundary condition set, including the ele- 
ments and nodes affected. 

5. Boundary solution (on an element basis), including displacements and 
tractions at each node of each element. The resultant load on each 
element and on the entire GMR is calculated and printed. 

6. Displacement, stress and strain on a nodal basis, at all surface 
nodes, for each GMR, 

7. Displacements at all volume cell nodes. 

8. Stresses at all volume cell nodes. 

9. Strains at all volume cell nodes. 
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Items 5 through 9 are printed at each user-requested time point. In the case 
of a plastic analysis, an extra print of the elastic solution at first yield 
is also provided. 

Sample output pages are included in Section 3.4.7, following the input de- 
scription. 

3.4.5 Program Validation/Verification 
3.4. 5.1 Validation of Elastic Capabilities 

The newly developed three-dimensional boundary element program (BEST) has a 
great variety of input and analysis options. More will be added in the future. 
It is clearly impossible to design and execute a test case examining every 
possible combination of these options. BEST has, however, been subjected to an 
extensive validation effort. 


A simplified tabulation of the elastic validation cases which have been used 
is presented in Figure 3.4-6. Major program capabilities are shown along the 
diagonal boundary. The box at the intersection of a row and column represents 
a possible test case combining two capabilities. The presence of an L, Q or M 
in such a box signifies that a test case has been successfully run using lin- 
ear, quadratic or mixed variation of displacement and traction. In many cases 
a single test problem involves several capabilities. This is indicated by mul- 
tiple entries in Figure 3.4-6. In most cases the test problems have been run 
repeatedly during the BEST development effort, in order to verify the con- 
tinued correct operation of existing capabilities after the addition of code 
enhancements. 
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Figure 3.4-6 BEST - Elastic Validation Matrix 
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A small subset of the validation cases, testing particularly important or 
basic capabilities, is discussed below. Except as otherwise indicated, all 
stresses are surface stresses, calculated using the displacements and trac- 
tions on the boundary. Unless otherwise noted, displacements are given in 
inches and stresses in Ib/in^. in many cases both full and hidden line plots 
of the meshes are shown. 

1. Cube in Tension (Figure 3.4-7) 


A cube was modeled using five quadrilateral and two triangular elements. 
The cube was subjected to a uniform displacement. Using a sufficiently 
accurate numerical integration, it is possible to achieve five place 
agreement with the exact solution. 
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Figure 3.4-7 Cube in Tension 



2. Pressurized Thick Cylinder (Figure 3.4-8) 

A 22.5 degree sector of a thick cylinder was modeled using six quadrilat- 
eral surface patches. In this problem two of the elements (the cylinder 
ID and OD) possess curvature. The problem was solved using both linear 
and quadratic variation of displacement and traction. The displacements 
from both analyses showed good agreement with analytical results. Quad- 
ratic variation was, however, required to obtain good stress results at 
the surface nodes. The use of more surface patches in the radial direc- 
tion would further improve the agreement between the calculated and anal- 
ytical values of the radial stress. 

This problem also tests the use of local coordinate systems, since the 
load is defined as a pressure (traction normal to the surface) on the ID 
of the cylinder, and rollers (zero normal deflection) are imposed on the 
0 and 22,5 degree planes. 
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Figure 3.4-8 Pressurized Thick Cylinder 



3. Rotating Thick Cylinder (Figure 3.4-9) 


The geometry and mesh for this problem are identical to (2), above. Note 
that for this problem the use of quadratic variation is required to ob- 
tain acceptable results for either displacement or stress. 
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Figure 3.4-9 Rotating Thick Cylinder 


4. Multi -GMR Cylinder (Figure 3.4-10) 

In this problem a 22.5 degree sector of a thick cylinder was again anal- 
yzed. In this case, however, several added features are present. First, 
the mesh is larger (twelve elements rather than six). Second, the struc- 
ture is broken into two GMRs along the 11.25 degree plane. Third, the 
boundary conditions on the 0 and 22.5 degree planes are cyclic symmetry 
between these two surfaces, rather than the roller boundary conditions 
used in (2) and (3), In this case only the results for quadratic varia- 
tion of displacements and tractions are shown. Once again, good agreement 
between the calculated and analytical results was obtained. 
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Figure 3.4-10 Thick Cylinder - Multi -GMR 


5. Half-Space Loaded with a Circular Punch (Figure 3.4-11) 

This problem was designed to test both the infinite elements incorporated 
in BEST and the ability to calculate stress and displacement at points in 
the interior of a region. 

The problem is the classical one of a half-space loaded with a uniform 
pressure over a circular region. In the figure the dashed outer lines of 
the mesh are boundaries at infinity. Use of the mesh shown, consisting of 
twelve finite quadratic elements and four infinite elements, gave good 
agreement between calculated and analytical results for the displacement 
and stress under the center of the punch. 
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Figure 3.4-11 Half -Space Loaded With a Circular Punch 


6. Long Beam on an Elastic Half-Space (Figure 3.4-12) 

In this problem a long beam is attached to an elastic half-space and 
loaded in tension. This problem is of interest primarily as a starting 
point for later dynamic/transient analyses. Note that the correct solu- 
tion for the tip deflection of the beam is essentially the sum of the ex- 
tension of the beam in simple tension and the displacement of the half- 
space under a patch load. 

It is of particular interest to note that it was never possible to obtain 
an acceptable solution to this problem when it was run as a single re- 
gion, using either linear or quadratic variation. This is due to the fact 
that, when considered as a single region, much of the beam is effectively 
located at infinity. To obtain accurate results would require a much more 
extensive mesh on the surface of the half-space, losing all the advantage 
gained by the use of infinite elements. Assigning the beam and the half- 
space to separate GMRs eliminates this problem, leading to reasonable re- 
sults (-5% error) for linear variation and very good results (less than 
Vio error) for quadratic variation. 
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Figure 3.4-12 Beam on an Elastic Half-Space 


3.4. 5.2 Validation of Inelastic Algorithms 

The main portions of the code which are unique to inelastic analysis are the 
constitutive modules and the routines for calculating volume integrals of var- 
ious kernel function/inelastic strain (or stress) products over twenty node 
isoparametric cells. All of the other basic functions of the program are iden- 
tical with the elastic version and do not, therefore, require separate valida- 
tion. Several simple problems have been run to verify the unique features of 
the inelastic code. 

1. Cube under Uniform Thermal Expansion 

An unconstrained unit cube was subjected to a uniform thermal expansion. 
Since the cube is unconstrained, the final stress should be zero and the 
final direct strains should be those due to the temperature change. To 
obtain these results using the inelastic code requires that all of the 
inelastic volume integrals be properly implemented. The problem was run 
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twice. In the both cases the surface was modeled using six (quadratic) 
elements. In the first analysis, a single volume cell was used. In the 
second, the cube was divided into four cells, requiring calculations at 
interior, as well as surface, points. 

The results were excellent for both analyses. 

2. Cube in Simple Tension (Figure 3.4-13) 

A cube in simple tension, represented by six surface elements and one 20 
noded cell was loaded in simple tension. The cube was modeled as an elas- 
tic linearly strain hardening plastic body. The analysis was carried out 
up to the fourth increment beyond the yield point. Excellent agreement 
with analytical results were obtained. The same problem was also analyzed 
with 12 surface elements and 4 volume cells (in order to test the volume 
integration routines for neighboring cells); the results were found to be 
essentially similar for both cases. 


y 




Figure 3.4-13 Elastic-Plastic Cube in Tension 
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3. Pressurized Thick Cylinder (Figures 3.4-14 through 3.4-17) 

An internally pressurized thick cylinder (inner radius a=l , outer radius 
b=2.0) was modeled using six surface elements and one cell. Axial deflec- 
tion on the front and back faces was constrained to simulate plane 
strain. Three-dimensional load-deflection results are compared to plane 
strain results in Figure 3.4-14. Even with this simple model, good ac- 
curacy was achieved over a significant range of nonlinearity when com- 
pared with the two-dimensional model. 
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Figure 3.4-14 Load Displacement Response for One Cell 
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Figure 3.4-15 shov/s a similar comparison with the two-dimensional results 
for a discretization using ten surface elements and two volume cells. 
Since the two-dimensional program is probably the most accurate boundary 
element method program analysis to date, it serves to establish the cor- 
rectness of the three-dimensional analysis. 
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Figure 3.4-15 Load Displacement Response for Two Cells 
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Figure 3.4-16 shows the convergence of the radial displacement solution 
when the plastic front has reached the middle of the cylinder (r=1.5a) 
for the three-dimensional solution, where it can be seen that even a two- 
cell (10 boundary elements) representation provides a satisfactory solu- 
tion. Similar convergence studies on the overall load-displacement re- 
sponse are shown in Figure 3.4-17, where the two-cell results are almost 
indistinguishable from the exact analytical solution. 



Figure 3.4-16 Variation of Radial Displacement When the Plastic Front Is At 
r«1.5a for a Three-Dimensional Analysis 
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Figure 3.4-17 Convergence of the Load Displacement Response for the Three- 
Dimensional Analysis of a Thick Cylinder 


4. Perforated Plate in Axial Tension (Figures 3.4-18 through 3.4-20) 

Figure 3.4-18 shows the discretization for a perforated plate under 
axial loading (the X-direction). Thirty boundary elements (8 noded) 
and six cells (20 noded) have been used to define the problem. It is 
important to note that the cells are defined only in the high stress 
concentration region where yielding is likely to occur. 

Figure 3.4-19 shows the overall load-displacement behavior obtained 
from the three-dimensional analysis compared with plane stress bound- 
ary element method results. Figure 3.4-20 shows the same comparison 
for the longitudinal stress distribution at the root of the plate. 
The results of the analyses agree very well with each other, indi- 
cating that the numerical implementation of nonlinear analysis within 
BEST is at least as accurate as the existing two-dimensional boundary 
element program. 
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Figure 3.4-18 Perforated Plate in Tension 
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3. 4. 5. 3 Validation of Dynamic Analysis 

The dynarnic/transient capabilities of the present program are implemented 
using transform techniques. The problem is recast in the Laplace or Fourier 
transform domain. In the case of a time harmonic loading the boundary element 
algorithm is exactly analagous to that for the static elastic case, except 
carried out in complex arithmetic. In the case of a more complex loading func- 
tion (either a discrete sum of harmonic terms or an aperiodic loading), an ap- 
propriate transform of the boundary condition is taken, and the system equa- 
tions are created and solved for a set of values of the transform parameter. 
The time domain solution must then be numerically reconstructed from the 
transform solution. 

Test cases have been developed to test the time harmonic solution, the numer- 
ical transform inversion and the ability to solve problems with nonharmonic 
loading. 

1. Cantilever Subject to Harmonic End Shear (Figure 3.4-21) 

A long cantilever was modeled using a total of eighteen quadratic surface 
patches. A time harmonic shear load at a frequency of 314 radians/second 
was applied to the free end of the beam. The excellent agreement between 
the calculated three-dimensional response and analytical results (based on 
beam theory) for the envelope of the vibrating beam is shown in the ref- 
erenced figure. 




Figure 3.4-21 Cantilever Subject to Harmonic End Shear 
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2. Cantilever Subject to Harmonic Transverse Load (Figure 3.4-22) 

The same model discussed in (1) was subjected to a time harmonic patch 
load. The agreement between the three-dimensional calculations and beam 
theory was, once again, excellent. 



Figure 3.4-22 Cantilever Subject to Harmonic Transverse Load 
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3. Test Problem for a Time Dependent Analysis (Figures 3.4-23 and 3.4-24) 

In order to test the numerical accuracy of the inversion of the 
transform domain solution, a problem of cantilever subjected to a 
time dependent loading at the free end was analyzed. 

Figure 3.4-23 shows the surface discretization used for the problem 
as well as the calculated displacement response. The end displacement 
agrees very well with the exact analytical solution. The difference 
is mainly due to the way the applied loading was represented by the 
direct transform algorithm for a piece-wise linear approximation, as 
shown in Figure 3.4-24. 
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Figure 3.4-23 Transient Analysis of a Cantilever Subjected to a Harmonic 
Axial Loading 
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Figure 3.4-24 Cantilever Loading - Actual vs. Simulated 


3. 4. 5. 4 Notched Specimen Verification 

Verification of this initial version of BEST is being carried out using test 
data and previous analytical results for notched specimens. The work done to 
date relates to specimen loading within the elastic range. The inelastic and 
creep analysis of one of these specimens is now being performed. 

1. C-Notch Low Cycle Fatigue Specimen (Figure 3.4-25) 

The C-notch low cycle fatigue specimen, as shown in the referenced figure, 
is designed to place a large volume of material in a plane strain, high 
stress condition. The specimen has, since its design several years ago, 
been subjected to a variety of elastic and inelastic analysis as well as 
to strain gage testing for specimen calibration. 
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Figure 3.4-25 Cross Section of C-Notch Low Cycle Fatigue Specimen 


In the present verification program the specimen was analyzed using BEST. 
The portion of the specimen which was analyzed is indicated in Figure 
3.4-25. The mesh used is shown, in both full and hidden line views, in 
Figure 3.4-26. The analysis was carried out using both linear and quad- 
ratic variation of displacement and traction. 

Key stress results are summarized in the table below. The plane strain 
results were obtained from a variety of two-dimensional codes. The base- 
line results were obtained from a very detailed three-dimensional analy- 
sis. The NASTRAN results cited were obtained using a mesh of twenty node 
isoparametric elements. The surface mesh refinement in the NASTRAN analy- 
sis was approximately equivalent to that in the BEST analysis. 
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Figure 3.4-26 BEST Model for Analysis of C-Notch Low Cycle Fatigue Specimen 


It is clear from the results above that BEST (using quadratic variation) 
is equivalent in accuracy to the previous baseline solution and to the 
plane strain results, and is superior to three-dimensional finite element 
analysis for an equivalent mesh. 

It is also clear that the use of linear variation for the full model does 
not provide sufficiently accurate results. The linear results could be 
improved by mesh refinement, but the use of quadratic variation over the 
same mesh is more efficient both in input preparation and analysis cost. 

The results of the quadratic BEST analysis are compared to strain gage 
data in Figure 3.4-27. The gages were located on the free surface of the 
specimen along the line shown in Figure 3.4-26. Agreement between the 
three-dimensional calculations and the data is generally good. 
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Figure 3.4-27 Comparison of Strain Gage Data With BEST Results for C-Notch 
Specimen 
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2. Benchmark Notch Specimen (Figure 3.4-28) 

The benchmark notch specimen is a double edge notch specimen developed by 
General Electric/Louisiana State University (GE/LSU) under NASA-Lewis 
Contract NAS3-22522 (Reference 83). A significant volume of well docu- 
mented data was provided as part of the referenced contract. These data 
have been used to verify the elastic capabilities of BEST. Verification 
of the inelastic capabilities of BEST, using this same data base, is now 
underway. 

The specimen geometry is defined in Figure 3.4-28. Stress analysis was 
carried out for the gage section only, a procedure already known to be 
satisfactory (Reference 83). Three different models were used, all ideal- 
izing one-quarter of the specimen gage section. Detailed comparison of 
results was carried out among the BEST analyses, GE/LSU finite element 
results and Pratt & Whitney finite element results obtained using both 
MARC and MARC-HOST. While these comparisons are not discussed here, it 
should be noted that, with sufficient mesh refinement, equivalent results 
were obtained with all analysis tools. 

The discussion in this report is directed at the comparison of BEST re- 
sults with the GE/LSU strain gage data. The major characteristics of all 
of the BEST analyses are summarized in the table below. The maximum per- 
ipheral strain in the notch (at the free surface) is given for each anal- 
ysis. All analysis methods and the test data show that this value should 
be between 1700 and 1800 (microstrain). 


Model 

Elements 

Variation 

Equations 

GMRs 

CPU Time 
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1 inear 
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64 
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50 

quadratic 
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mi xed 
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1 i near 
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10 
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10 

quadratic 

96 
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28 

1605 

It is 

clear from 

the table 

that models 

2, 4 and 

5 all yield 

results 


accuracy entirely comparable with the strain gage data. The variation in 
peak strain among these three is within +1.5%. It is also clear that the 
most cost effective analysis is that whTch combines substructuring with 
mixed linear and quadratic variation. As was the case for the C-notch 
specimen, fully linear analysis cannot give acceptably accurate results 
without unacceptable mesh refinement. 
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Figure 3.4-28 Double Edge Notch Specimen Used in Contract NAS3-22522 


The mixed model is shown in Figure 3.4-29, in both full and hidden line 
views. In this model eight of the twenty-two elements (those in and near 
the notch) were quadratic, while the remaining fourteen were linear. The 
visible elements in the hidden line plot are identified as linear or 
quadratic. The fully linear and quadratic models (4 and 6) utilized an 
identical surface discretization to that shown. The results of these 
three analyses (4, 5 and 6) are plotted with the GE/LSU strain gage re- 
sults in Figure 3.4-30. Both the fully quadratic and mixed analyses are 
in excellent agreement with the strain gage data. The difference betv^/een 
the two analyses is far less than the normal scatter in strain gage data. 
The fully linear analysis, however, does not give either an accurate peak 
strain or a correct representation of the strain distribution near the 
notch. 

Also shown in Figure 3.4-29 is a pattern of eight volume cells. These 
cells are used for the representation of inelastic strain in the nonlin- 
ear analysis of the specimen under various loadings. The cell pattern was 
designed to include the portion of the specimen gage section which exper- 
iences inelastic response during the tests described in Reference 83. 
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Figure 3.4-29 Optimized Model for Analysis of Benchmark Notch Specimen 
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Figure 3.4-30 Comparison of BEST Results With Strain Gage Data 
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3.4. 5. 5 Hot Section Component Analysis 

In order to evaluate the capabilities of the present version of BEST for the 
analysis of real components, an analysis of a commercial cooled turbine blade 
geometry was carried out. It is expected that the use of this analysis as a 
benchmark problem will be continued throughout the life of the Inelastic Meth- 
ods program. This report discusses preliminary results for the elastic analy- 
sis of the blade. 

The blade analyzed is a cooled high turbine blade presently in service. It is 
subject to mechanical loads (primarily centrifugal) and thermal loads. Of par- 
ticular interest for this blade is the location and magnitude of the peak 
stress under the platform. 

A BEST model was built for this problem. Both full and hidden line views of 
the model are shown in Figure 3.4-31. The model consists of five GMRs. The 
interfaces between GMRs are generally perpendicular to the radial direction. 
The characteristics of the model are summarized below. 


GMR 

El ements 

Linear 

Source Points 

Quadratic 
Source Points 

1 

60 

61 

180 

2 

86 

85 

256 

3 

107 

98 

303 
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106 

96 

298 

5 

80 

76 

232 


The system equations for a fully linear analysis contain 1248 equations. 

To date only a fully linear analysis for centrifugal loads has been carried 
out. Analysis of the results is still in a preliminary stage. The total cen- 
trifugal load at various spanwise stations on the blade has been compared with 
the design calculations for the blade (Figure 3.4-32). The agreement between 
the two totally independent calculations is excellent. The total centrifugal 
load for the blade is within 1 % of the design calculation. The larger local 
error near the blade platform is believed to be due to the fact that the de- 
sign calculation models the platform as a discontinuous addition of mass. 

Extremely preliminary study of blade tip deflections, load distribution over 
the base of the blade neck and concentrated stresses indicates reasonable 
agreement with finite element results. Contour plots of principal stress con- 
tours have verified that the peak stress occurs in the correct location. It is 
clear, however, from the results of the verification problems that at least 
some use of quadratic variation will be required to achieve correct definition 
of critical stresses and strains, 

A fully quadratic analysis of this problem would lead to a system of over 3800 
equations. Analysis time, on the IBM 3081, would be approximately one hour, 
compared to about 16 minutes for the fully linear analysis. The results of the 
benchmark notch analysis clearly indicate that it should be possible to 
achieve acceptably accurate results at reasonable cost through the use of 
mixed variation. Present efforts are directed at identifying those areas of 
the model requiring quadratic variation and/or mesh modification. 
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Figure 3.4-31 BEST Model of Cooled Turbine Blade 


3.4-57 


CENTRIFUGAL LOAD (lb) 
THOUSANDS 



% SPAN 

DESIGN CALCULATION + BEST 


Figure 3.4-32 Centrifugal Load Distribution in a Cooled Turbine Blade 
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3.4.6 Boundary Element Stress Technology (BEST) Program U911 - Input De- 
scription 

The input to U911 is divided into six sections as follows: 

1. Case control (**CASE control) 

2. Material properties (**MATErial property) 

3. Generic modeling regions (**GMRegion) 

4. Interfaces (**INTErface) 

5. Boundary condition sets (**BCSEt) 

6. Body forces (**B0DY force). 

A detailed description of each of these sections is provided in the following 
paragraphs. The interface and body force sections are optional; the other sec- 
tions must be input at least once. 

Input quantities may be either alphanumeric or numeric (integer, floating 
point, E, or D format) as specified and may be up to 16 characters. Individual 
entries on a card (both keyword and input) must be separated by at least one 
blank space. Input for certain keywords (as noted) may be continued onto more 
than one card by repeating the keyword on the new card(s). 

Keywords may be input as shown; minimum input is the CAPITALIZED characters. 
Those keywords which are underscored must always be input. Keywords shown be- 
low are indented to indicate groups of cards to be input together. However, it 
is not necessary to indent the input in this manner. 

Consistent units must be used. Angles are in degrees, speed is in revolutions 
per minute, and frequency is in radians per sound. 

The current program limits include: 

20 time points 

10 generic modeling regions 
15 surfaces per region 

600 elements (416 quadratic elements in problems with body forces, 
300 elements in problems having interior points) 

11 infinite elements per region 
2500 nodes (560 nodes per region) 

1200 source points (600 source points in problems having interior 
points) 

302 source points per region in a local coordinate system 
99 interface and cyclic symmetry element pairs (total) 

350 interface and cyclic symmetry node pairs (total) 

20 boundary condition sets with cyclic symmetry 

60 boundary condition sets with springs 

maximum element number of 9999 

maximum node number of 9999 

maximum of 24 entries per input card. 
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3. 4. 6.1 Case Control Input 

The case control input section may be input only once. 



Type of 


Keyword 

Input 

Input 

**CASE control 



TTTre 

miTri fugal 

A1 phanumeri c 

Case title 

DYNAmic 

Numeric 

Frequency value, damping coefficient 

INHOmogeneous 

THERmal 



P LAS ti city 
RESTart 

A1 phanumeri c 

READ or WRITE 

TIMES 

Numeric 

Output time value(s) 

TRANsient 

Numeri c 

Number of time intervals, time 
increment, damping coefficient 

The analysis is 

assumed to be static 

, homogeneous, constant temperature. 

elastic, and time 

independent unless the appropriate optional keyword is in- 

put. The optional 

keywords need be included only if a particular option is to 

be turned on. 



The case title should have a maximum of 72 characters. 

Input on the TIMES 

card may be continued 

on more than one card. 

3. 4. 6. 2 Material 

Property Input 


The material property input section 

must be repeated for each separate 

material , 

Type of 


Keyword 

Input 

Input 

**MATErial property 


ID 

ISOTropic 

A1 phanumeri c 

Material name 

Temperature value(s) 
Young's modulus value(s) 

TEMPerature 

Numeric 

EMODulus 

Numeric 

ALPHa 

Numeri c 

Alpha value(s) 

DENSity 

Numeric 

Mass density value 
Poisson's ratio value 

POISson 

Numeric 
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Keyword 

Type of 
Input 

Input 

AN isotropic 

TEMPerature 

Numeri c 

Temperature value 

ALPHa 

Numeri c 

Alpha value 

CONStants 

Numeri c 

Elastic constant values 

DS 

INELastic 

MONOtonic 

CYCLic 

ISOTropic 

TIME 

Numeri c 

Time point identifier(s) for slow 

YIELd 

Numeric 

algorithm 

Proportional limit value 

CURVe 

Numeri c 

Stress value, plastic strain value 

TWO surface 
TIME 

Numeric 

Time point identifier(s) for slow 

YIELd 

Numeric 

al gori thm 

Inner proportional limit value. 

HARD 

Numeric 

outer proportional limit value 
Inner hardening parameter, outer 

WALKER 


hardening parameter 


If the material name is a material library designation, then the TEMPerature, 
EMODulus, ALPHa, DENSity, and POISson cards may be omitted. However, any of 
them may be included to override the material library data. 

The Young's modulus and alpha values must be input in the same order as the 
temperature values. 

A DS card (indicating a directionally solidified material) must be included in 
ANISotropic input. The constants to be input for a directionally solidified 
material are C]] , C-]3, C33, C44, and Cse- 

Either a MONOtonic card (indicating monotonic loading) or a CYCLic card (in- 
dicating cyclic loading) must be included in INELastic input. Either the 
ISOTropic model, ^ the TWO surface model, £r the WALKer must be used in 
INELastic input. 

The value of stress at zero plastic strain (i.e., proportional limit) need not 
be included in CURVe input. The stress and strain values input should start 
with the first nonzero value of plastic strain, A HARDening parameter is cal- 
culated and used if a single stress/strain pair is input. Otherwise, the mul- 
tipoint stress/strain algorithm is used. 

Input on the TEMPerature, EMODulus, ALPHa, CONStants, and YIELd cards may be 
continued on more than one card. 
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3. 4. 6. 3 Generic Modeling Region Input 

The generic modeling region section must be repeated for each region. 


Type of 

Keyword Input Input 


**GMR 


ID 

A1 phanumeric 

TRHOmogeneous 


THERmal 


PLAStici ty 


MAT 

A1 phanumeric 

TREFerence 

Numeric 

POINts 

Numeric 

SURFace 

A1 phanumeric 

TYPE 

Alphanumeric 

ELEMents 

Numeric 

TRANslate 

Numeri c 

REF 

Numeri c 

DIR 

Numeri c 

ROT 

Numeric 

NORMal 

Alphanumeric 

VOLUme 


POINts 

Numeric 

TYPE 

A1 phanumeric 

CELLS 

Numeric 


Region name 


Material name 

Reference temperature value 
Node number, coordinate values 
{x,y,z) 

Surface name, (reference surface 
name) 

LINI or QUAD 

(I), "element number, node numbers, 
(ref, node) 

Translation values (x,y,z) 

Axis of rotation reference point 
(x,y,z) 

Axis of rotation direction (x,y,z) 
Axis of rotation angle (degrees) 
Element number, + or - 

Point number, coordinate values 
(x,y,z) 

LINI or QUAD 

Cell number, node numbers (nodes 1 
to 10) 

Cell number, node numbers (nodes 11 
to 20) 


The SURFace input may be either of two forms: 

- a TYPE card and an ELEMents card to define element connectivity 

- a REF card, a DIR card, a ROT card, and/or a TRANslate card to define 
one surface with reference to another surface (translations are per- 
formed first, followed by rotations). 

The TYPE designation in SURFace input specifies the traction or displacement 
variation on the element. A surface may contain only one TYPE card. Therefore, 
if mixed variation is required in a region, two surfaces must be defined. 


3.4-62 



Elements must have either 6 (triangles) or 8 (quadrilaterals) nodes. Element 
numbering is consecutive around the boundary. 

Infinite elements are indicated by an I on the element card and may have 7 
(triangles) or 9 (quadrilaterals) nodes, where the extra node is the reference 
node. If the reference node is not input, it is assumed to be at the origin. 

The sign associated with the defining element on the NORMal card should be 
plus (+) if the element is numbered in a counterclockwise direction as seen 
from the outside of the model or minus (-) if it is numbered in a clockwise 
direction. Disjoint boundaries must have multiple element/sign pairs on the 
NORMal card. 

The points which are input in the VOLUme input (both points and cell nodes) 
are treated as "interior" points. These points may be either nodal points, 
other surface points, or true interior points. The TYPE designation in VOLUme 
input specifies cell source points (i.e., 8 corner nodes or all 20 nodes). 

Cells must have 20 nodes. Cell numbering is consecutive around the "front" 
face boundary, followed by the four midplane nodes, followed by the "back" 
face boundary. 

3. 4. 6. 4 Interface Input 

The interface section must be repeated for each interface. 


Keyword 

Type of 
Input 

Input 

**INTErface 

GMR 

Alphanumeric 

Region name of first region 

SURFace 

Alphanumeric 

Surface name in first region 

ELEMents 

Numeri c 

Element number(s) in first region 

POINts 

Numeric 

Node number (s) in first region 

GMR 

Al phanumeric 

Region name of second region 

SURFace 

A1 phanumeric 

Surface name in second region 

ELEMents 

Numeric 

Element number (s) in second region 

POINts 

Numeric 

Node number(s) in second region 

SLIDing 
The interface 

is assumed to have complete 

displacement compatibility unless a 


SLIDing card is input, in which case only normal displacement compatibility is 
assumed. 

The ELEMents card and/or the POINts card are included in SURFace input only to 
designate a subset of that surface. 

Input on the ELEMents and POINts cards may be continued on more than one card., 
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3.4. 6. 5 Boundary Condition Set Input 

The boundary condition set section must be repeated for each new boundary con- 
di ti on . 



Type of 


Keyword 

Input 

Input 

**BCSEt 

ID 

A1 phanumeric 

Boundary condition set name 

GMR 

VaTiiq 

Al phanumeric 

Region name 

1 riL. 

RELAtion 

SURFace 

Alphanumeric 

Surface name 

ETEMents 

Numeric 

Element number(s) 

POINts 

Numeric 

Node number (s) 

TIMES 

Numeric 

Input time value(s) 

LOCAl 

CYCLic 

GMR 

Al phanumeric 

Region name 

SURFace 

Al phanumeric 

Surface name 

ELEMents 

Numeric 

Element number (s) 

POINts 

Numeric 

Node number(s) 

ANGLe 

Numeric 

Axis of rotation angle (degrees) 

DIR 

Numeri c 

Axis of rotation direction (x,y,z) 

Displacement 

Numeric 

Component value 

SPLIst 

Numeric 

Source point value(s) or ALL or SAME 

T 

Numeric 

Time point identifier, displacement 



value(s) 

RIGId 

Numeri c 

Component value 

SPRIng 

Numeric 

Component value, spring value 

TRACtion 

SPLIst 

Numeric 

Source point value(s) or ALL or SAME 

T 

Numeric 

Time point identifier, traction 
value(s) 


The ELEMents card and/or the POINts card are included in SURFace input only to 
designate a subset of that surface. 

The TIMES card must be omitted in a boundary condition set which contains a 
RIGId card. If the value(s) on the TIMES card differ from those values in the 
case control input, the output is calculated by linear interpolation. In the 
case of time independence (i.e,, the TIMES card is omitted) the time point 
identifier on the T card must be 1 (one). 

The LOCAl card designates input in the outward normal direction. The component 
value on the CYCLic, Displacement, RIGId, SPRIng, or TRACtion card must be 1 
(one). Care must be taken not to mix global and local coordinate systems on a 
particular element. Care must also be taken not to input conflicting com- 
ponents on a particular node in a particular element. 
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The VALUe card should be included with the Displacement card, the RIGId card, 
or the TRACtion card. The RELAtion card should be included with the CYCLic 
card or the SPRIng card. 

Either CYCLic input, or Displacement input, or RIGId input, ^ SPRIng input, 
or TRACtion input must~5e included in a boundary condition set. This input set 
may be included up to three times (once for each component) in a boundary con- 
di ti on set. However, different boundary condition types may not be mixed in a 
boundary condition set. 

The cyclic symmetry direction (DIR card) defaults to the z-axis (0,0,1). 

The SPLIst card indicates the order in which the values are to be input on the 
T cards. The input may be in either of three forms: 

- nodal values 

- ALL to indicate that a single constant value is to be input 

- SAME to use the previous source point list within the current boundary 
condition set (this option may not be used for the first source point 
list in the current boundary condition set). 

Input on the ELEMents, POINts, SPLIst, and TIMES cards may be continued on 
more than one card. Input on the T card may be continued on more than one 
card, including the time point identifier on each card. 

3. 4. 6. 6 Body Force Input 

The body force section is optional. 



Type of 


Keyword 

Input 

Input 

**B0DY force 



CENTrifugal 



DIR 

Numeric 

Axis of rotation direction (x,y,z) 

PT 

Numeric 

Axis of rotation reference point 
(x,y,z) 

SPEEd 

Numeric 

Speed value(s) 

TIMES 

THERmal 

Numeric 

Time value(s) 

TIMES 

Numeric 

Time value(s) 

TEMPeratures 

Numeric 

Node number, temperature value (s) 
at time(s) 

The CENTrifugal 

card and the TIMES card must be included in the case control 

section in order 

to input centrifugal 

body forces. 


The direction (DIR card) defaults to the z-axis (0,0,1). The reference point 
(PT card) defaults to the origin (0,0,0). 

Input on the SPEEd and TIMES cards may be continued on more than one card. 
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3.4.7 Sample Output from BEST 
3. 4. 7.1 Input Echo 


INPUT ECHO )*»»» 


LINE INPUT 

1 »*CASE 

Z TITLE U911 TEST CASE - CUBE IN SIHPLE TENSION WITH PLASTICITY 

3 TIHES 1.00 2.00 3.00 4.00 5.00 

4 PLASTICITY 

5 »*MATE 

6 ID MATl 

7 TEMP 70.0000 

6 EMOD 100.000 

9 POIS .250000 

10 DENS 1.00000 

11 INELASTIC 

12 YIELD 100. 

13 TIMES 4 

14 CURVE 

15 150. 1. 

16 »»GMR 

17 ID GMRl 

18 PLASTICITY 

19 MAT MATl 

20 TREF 70. 

21 POINTS 


22 

1 


,0 


.0 


.0 


23 

2 


.0 


.0 


.500000 


24 

3 


.0 


.0 


1.00300 


25 

11 


.0 


.500000 


.0 


26 

13 


.0 


.500000 


1.00000 


27 

21 


.0 


1.00000 


.0 


28 

22 


.0 


1.00000 


.500000 


29 

23 


.0 


1.00000 


1.00000 


30 

101 


.500000 


.0 


.0 


31 

103 


.500000 


. 0 


1.00000 


32 

121 


.500000 


1.00000 


.0 


33 

123 


.500000 


1.00000 


1.00000 


34 

201 


1.00000 


.0 


.0 


35 

202 


1.00000 


.0 


.500000 


36 

203 


1.00000 


.0 


1.00000 


37 

211 


1.00000 


.500000 


.0 


38 

213 


1.00000 


.500000 


1.00000 


39 

221 


1.00000 


1.00000 


.0 


40 

222 


1.00000 


1.00000 


.500000 


41 

223 


1.00000 


1.00000 


1.00000 


42 

SURFACE SURFll 






43 

TYPE LINI 







44 

ELEMENTS 







45 


1 

1 

2 

3 13 

23 

22 21 

11 

46 


2 

201 202 

203 213 

223 

222 221 

211 

47 


3 

1 

2 

3 103 

203 

202 201 

101 

46 


4 

21 

22 

23 123 

223 

222 221 

121 

49 


5 

1 

11 

21 121 

221 

211 201 

101 

50 


6 

3 

13 

23 123 

223 

213 203 

103 

51 

NORMAL 


1 + 







52 VOLUME 

53 TYPE QUAD 

54 CELLS 
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3, 4. 7. 2 Case Control Summary 


CASE CONTROL INPUT 

JOB TITLE U911 TEST CASE - CUBE IN SIMPLE TENSION WITH PLASTICITY 


TIMES FOR SOLUTION: 

1.00 2.00 

3.00 4.00 


5.00 

PLASTICITY FLAG 

1 

THERMAL FLAG 

0 


CENTRIFUGAL LOAD FLAG 0 

TRANSIENT FLAG 

0 


DYNAMIC FLAG 0 


INTEGRATION EPSILON 

0.0010 

INTEGRATION GRADING 

FACTOR 1.4142 

RESTART FLAG 

0 



»»»» MATERIAL INPUT *»»» 

MATERIAL NAME MATI 

INELASTIC 

ISOTROPIC 

YIELD STRESS: 100.00 

STRESS STRAIN 

O.lOOOOE+03 0.0 

0.I5000E+03 O.IOOOOE+OI 

DENSITY: 1.0000 POISSONS RATIO: 0.E500 

TEMP ALPHA E 

0,70000E+Oa 0.0 O.IOCOOE+03 


3. 4, 7. 3 Generic Modeling Region (GMR) Definition 


«»»» gmr input 


REGION 

1 










MANE 6MR1 

MATERIAL MATI 


REFERENCE TEMPERATURE 

70.00 


NODES 

2Ct 

ELEMENTS 6 


SURFACES 

1 



SOURCE 

POINTS 

8 CELLS 

1 


INFINITE ELEMENTS 0 



COORDINATE LIST 










NODE 

X 

Y 


Z 


NODE 

>: 

Y 

Z 

1 

0.0 

0.0 


0.0 




0.0 

C.O 

0.5000 

3 

0.0 

0.0 


1.0000 



11 

0.0 

C.SC'OC 

0.0 

13 

0.0 

o.sooo 


1.0000 



21 

0.0 

1 .cooo 

0.0 

ZZ 

0.0 

1.0000 


0.5000 



23 

0.0 

1 .cc;o 

l.OOOO 

101 

0.5000 

0.0 


0.0 



103 

0.50CO 

0.0 

l.OOOO 

121 

0.5QQ0 

1.0000 


0.0 



123 

0.5000 

l.COCO 

1.0000 

201 

1.0000 

0.0 


0.0 



202 

1.0000 

C.O 

0.5000 

203 

1.0000 

0.0 


1.0000 



211 

1.0000 

o.ccoc 

D.O 

213 

l.OOOO 

0.5000 


1.0000 



221 

1.0000 

l.OOOO 

0.0 

222 

1.0000 

1.0000 


0.5000 



223 

1.0000 

l.OCCO 

1.0000 

SURFACE 

SURFU 

LINEAR 

VARIATION 








1 

1 

2 

3 

13 

23 

22 

21 

a 



2 

201 

211 

221 

222 

223 

213 

203 

2C2 



3 

1 

101 

201 

202 

203 

103 

3 

2 




21 

22 

23 

123 

223 

222 

221 

121 



5 

1 

11 

21 

121 

221 

211 

201 

101 



6 

3 

103 

203 

213 

223 

123 

23 

13 



SOURCE POINT LIST 


3 21 

23 

201 

203 

221 

223 













CELL INPUT 
3 2 

1 

11 

21 

22 

23 

13 

103 101 

121 

123 

203 

202 

201 

211 

221 

222 

223 

213 
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3. 4. 7. 4 Boundary Condition Definition 


BOUNDARY CONDITION INPUT »*»* 

BOUNDARY CONDITION SET NAME TRACll TYPE: VALUE 

GMR GMRl SURFACE SURF 11 

ELEMENT LIST 
Z 

SOURCE POINT LIST 

201 203 223 221 

LOCAL (NORMAL) COORDINATE SYSTEM 


TIME VALUES: 

1.0000 2.0000 3.0000 4.0000 5.0000 


COMPONENT 1 TRACTION INPUT 


TIME = 1.00 DATA VALUES: 

O.lOOOOE+03 O.lOOOOE+03 O.lOOOOE+03 O.lOOOOE+03 

TIME = 2.00 DATA VALUES: 

0.10500E+03 0.10500E+03 0.10500E+03 0.10500E+03 


TIME = 3.00 DATA VALUES: 

0.12000E+03 0.12000E+03 0.12000E+03 0.12000E+03 

TIME = 4,00 DATA VALUES: 

0.11250E+03 0.112SOE+03 0.11250E+03 0.11250E+03 

TIME 5 5.00 DATA VALUES: 

0.67S00E+O2 0.67500E+02 0.67500E+02 0.67500E+02 


3.4, 7, 5 Boundary Solution (Element Basis) 


JOB TITLE: U9U TEST CASE - CUBE IN SIHPLE TENSION WITH PLASTICITY 

BOUNDARY SOLUTION FOR TINE « t.OO 


ELEMEKT 

NODE NO. 

X OISPLACEMEWT 

Y OISPUCEHEKT 

Z DISPLACEMENT 

X TRACTION 

Y TRACTION 

Z TRACTION 

1 

1 

0.0 

0.0 

0.0 

'0.200000403 

0.0 

0.0 

1 

3 

0.0 

0.0 

-0.25001D>00 

>0.100000^03 

0.0 

0.0 

1 

23 

0.0 

•0.25002D«00 

•0.2S001D«00 

>0.10000D«03 

0.0 

0.0 

1 

21 

0.0 

«0. 250020^00 

0.0 

>0.100000^03 

0.0 

0.0 


2 

201 

0.10000D«01 

0.0 

0.0 

0.100000403 

0.0 

0.0 

2 

221 

0.100000«01 

>0.250010409 

0.0 

0.100000403 

0.0 

0.0 

2 

223 

O.lOOOOO^Ol 

>0. 250020409 

>0.250020400 

0.100000403 

0.0 

0.0 

2 

203 

0.100000401 

0.0 

>0.250020400 

0.100000403 

0.0 

0.0 


3 

1 

0.0 

0.0 

0.0 

0.0 

>0.429690-05 

0.0 

3 

201 

0.100000401 

0.0 

0.0 

0.0 

0.146930-04 

0.0 

3 

203 

0.10000D401 

0.0 

>0.250020400 

0.0 

-0.117610-04 

0.0 

3 

3 

0.0 

0.0 

>0.250010400 

0.0 

0.632120-06 

0.0 


4 

21 

0.0 

-0.250020400 

0.0 

0.0 

0.0 

0.0 

4 

23 

0.0 

>0.250020400 

>0.250010400 

0.0 

0.0 

0.0 

4 

223 

0.100000401 

>0.250020400 

-0.250020400 

0,0 

0.0 

0.0 

4 

221 

0.100000401 

>0.250010400 

0.0 

0.0 

0.0 

0.0 


5 

1 

0.0 

0.0 

0.0 

0.0 

0.0 

-0.205750-04 

5 

21 

0.0 

-0.25002D400 

0.0 

0.0 

0.0 

-0,276970-04 

5 

221 

0.100000401 

-0.250010400 

0.0 

0.0 

0.0 

0.216390-04 

5 

201 

0.100000401 

0.0 

0.0 

0.0 

0.0 

0.362410-04 


6 

3 

0.0 

0.0 

-0.250010400 

0.0 

0.0 

0.0 

6 

203 

0.100000401 

0.0 

>0.250020400 

0.0 

0.0 

0.0 

6 

223 

0.100000401 

-0.250020400 

>0.250020400 

0.0 

0.0 

0.0 

6 

23 

0.0 

-0.250020400 

>0.250010400 

0.0 

0.0 

0.0 
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3. 4. 7. 6 Boundary Solution (Nodal Basis) 


JOB TITLE: U911 TEST CASE - CUBE IN SIMPLE TENSION WITH PLASTICITY 

NODAL OUTPUT AT TIME = 1.00 

NODE DISPLACEMENT STRESS STRAIN 

X/Y/2 XX/YY/22 XY/X2/Y2 XX/YY/22 XY/X2/Y2 


1 

0.0 

0.0 

0.0 

0.10000E'f03 

-0.11716E-02 

-0.11579E-02 

0.0 

0.0 

0.0 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.0 

0.0 

0.0 

z 

0.0 

0.0 

-0.12501E+00 

O.lOOOOE-^03 

-0.10033E-02 

-0.17395E-02 

0.0 

-0.23842E-05 

-0.23842E-05 

O.lOOOOE+01 

-0.25001E+00 

-0.2S001E+00 

0.0 

-0.29802E-07 

-0.29802E-07 

3 

0.0 

0.0 

-0.25001E+00 

0.10000E403 

-0.11691E-02 

-0.11698E-02 

0.0 

-0.15895E-05 

0.0 

O.lOOOOE+01 

-0.2500IE+00 

-0.25001E+00 

0.0 

-0.19868E-07 

0.0 

11 

0.0 

-0.12S01E+00 

0.0 

0.10000E-f03 

-0.17748E-02 

-0.99501E-03 

0.23842E-05 

0.0 

0.0 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.29802E-07 

0.0 

0.0 

13 

0.0 

-0.12501E+00 

-0.25001E+00 

O.lOOOOE'i'03 

-0.17519E-02 

-0.99945E-03 

0.0 

0.0 

0.23842E-05 

O.lOOOOE+01 

-0.25001E+C0 

-0.2S001E+00 

0.0 

0.0 

0.29S02E-07 

21 

0.0 

-0.25002E+00 

0.0 

0.10000E+03 
-0.11924E-02 
-0. 11606E-02 

0.31789E-05 

0.0 

0.1B89SE-05 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.39736E-07 

0.0 

0.19668E-07 

22 

0.0 

-0.25002E+00 

-0.12501E+00 

0.10000E-»03 

-0.10028E-02 

-0.17242E-02 

0.0 

-0.47684E-05 

0.23842E-05 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.0 

-0.59605E-07 

0.29802E-07 

23 

0.0 

-0.25002E+00 

-0.25001E+00 

O.lOOOOE+03 

-0.11628E-02 

-0.11546E-02 

-0.15895E-05 

-0.63578E-05 

0.31789E-05 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

-0.19866E-07 

-0.79473E-07 

0.39736E-07 

101 

O.SOOOOE-fOO 

0.0 

0.0 

O.lOOOOE+03 

-0.75028E-03 

-0.75923E-03 

-0.28610E-04 

0.0 

0.0 

O.lOOOOE+01 

-0.2S001E+00 

-0.25001E+00 

-0.35763E-06 

0.0 

0.0 

103 

0.50000E+00 

0.0 

-0.25002E+00 

0. 100001: -»03 
-0.744901--03 
-0.74768I--03 

-0.95367E-05 

-0.47684E-05 

0.0 

O.IOOOOE+Ol 

-0.25001E+00 

-0.25001E+00 

-0.11921E-06 

-0.59605E-07 

0.0 

121 

0.50000E+00 

-0.25001E+00 

0.0 

0.100001:^03 

-0.75531l;-03 

-0.76147I--03 

0.11921E-04 

0.19073E-04 

0.0 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.14901E-06 

0.23842E-06 

0.0 

123 

0.50000E+00 

-0.25002E+00 

-0.25002E+00 

0.100001;+03 

-0.74005I;-03 

-0.74005E-03 

0.26226E-04 

0.11921E-04 

0.0 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

0.32783E-06 

0.14901E-06 

0.0 

201 

O.lOOOOE+01 

0.0 

0.0 

O.lOOOOE+03 

-0.11649i;-02 

-0.11921E-02 

-0.50863E-04 

0.0 

0.0 

O.lOOOOE+01 

-0.25001E+00 

-0.25001E+00 

-0.63578E-06 

0.0 

0.0 
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3. 4. 7. 7 Cell Node Displacements 


JOB TITLE: U911 TEST CASE - CUBE IN SIMPLE TENSION WITH PLASTICITY 

INTERIOR DISPLACEMENT AT TIME =1.00 


NODE 

X DISPLACEMENT 

Y DISPLACEMENT 

Z OISPLACEMENT 

3 

0.0 

0.0 

-0.2500150+00 

2 

0.2900250-07 

-0.7951030-08 

-0.1250070+00 

1 

0.0 

0.0 

0.0 

n 

0.29802SD-07 

-0. 1250080+00 

-0.7951020-08 

El 

0.0 

-0.2500150+00 

0.0 

22 

0.298025D-07 

-0.2500150+00 

-0.1250070+00 

23 

0.0 

-0.2500150+00 

“0.2500150+00 

13 

0.2980250-07 

-0.1250080+00 

-0.2500150+00 

103 

0.5000020+00 

-0.7951030-08 

-0.2500150+00 

101 

0. 5000020+00 

-0.7951030-08 

-0.7951030-08 

121 

0.5000020+00 

-0.2500150+00 

-0.7951030-08 

123 

0.5000020+00 

-0.2500150+00 

-0.2500150+00 

203 

0.1000000+01 

0.0 

-0.2500150+00 

202 

0.1000000+01 

-0.7951030-08 

-0.1250080+00 

201 

0.1000000+01 

0.0 

0.0 

211 

0.1000000+01 

-0.1250060+00 

-0.7951090-08 

221 

0.1000000+01 

-0.2500150+00 

0.0 

222 

0.1000000+01 

-0.2500150+00 

-0.1250080+00 

223 

0.1000000+01 

-0.250015D+Q0 

-0.2500150+00 

213 

0.1000000+01 

-0.1250080+00 

-0.2500150+00 


3.4, 7.8 Cell Node Stresses 


JOB TITLE: U911 TEST CASE - CUBE IN SIMPLE TENSION WITH PLASTICITY 

INTERIOR STRESS AT TIME a I. 00 


NODE 

SIGMAXX 

SIGMAYY 

SIGMAZZ 

TAUXY 

TAUXZ 

TAUYZ 

3 

0.100000D<»03 

-0.U7113O-0Z 

-0.116Z72D-02 

0.0 

-0. 36984 7D-05 

0.1861500-07 

Z 

0.100000D403 

•0.101165D»0Z 

-0.1746280-02 

0.0 

-0.2773850-05 

0.1396120-07 

1 

0.1000000«03 

-0.117198D-0Z 

-0.1 162930-02 

0.0 

0.0 

0.0 

11 

0.1000000^03 

•0.175666O~0Z 

-0.9981750-03 

0.7072690-06 

0.0 

-0.4559900-06 

Z1 

0.1000000«03 

-0.U7Z93O-0Z 

-0.1161260-02 

0.9430530-06 

0.0 

-0. 6079870-06 

zz 

0.100000D«03 

>0.101Z690>0Z 

-0.1746740-02 

0.0 

-0. 2218120-05 

-0.8980190-06 

Z3 

0.1000000^03 

-0.1170980-OZ 

-0.1162630-02 

-0.6126240-06 

-0.2957490-05 

-0.5893720-06 

13 

0.100000003 

-0.17533ZO-0Z 

-0.1006150-02 

-0.6096160-06 

0.0 

-0.4280680-06 

103 

0. 1000000*03 

-0.745Z6ZO-03 

-0.7458010-03 

0.1470830-06 

-0.1792190-06 

0.0 

101 

0.1000000*03 

-0.7531630>03 

-0.7586700-03 

-0.3913120-06 

0.5366490-05 

0.0 

IZl 

0.1000000*03 

-0.7<»9Z9ZD-03 

-0.7513860-03 

0.1023270-05 

0.5906880-05 

0.0 

1Z3 

0.1000000*03 

-0. 74 76<t 80-03 

-0.7449760-03 

-0.107215D-05 

0.1470640-05 

0.0 

Z03 

0.1000000*03 

-0.117Z99D-0Z 

-0.1160600-02 

0. 1961110-06 

0.3459510-05 

0,7595940-06 

ZOZ 

0.1000000*03 

-O.1O19Z0O-OZ 

-0.1778620-02 

0.0 

0.7963120-05 

0,5696950-06 

ZOl 

0.1000000*03 

-O.U77990-OZ 

-0.1199010-02 

-0. 5217500-06 

0. 7157960-05 

0,0 

Zll 

0.1000000*03 

-0.1761350-OZ 

-0.1037160-02 

-0. 7533500-07 

0.0 

-0.1772900-05 

ZZl 

0.1000000*03 

-0. 1171340-02 

-0.1166660-02 

0.4213030-06 

0.7875840-05 

-0.2363860-05 

zzz 

0.1000000*03 

-0.101771D-OZ 

-0.1772660-02 

0.0 

0.9595640-05 

-0,2976100-05 

ZZ3 

0.1000000*03 

-0. 1178270-02 

-0.1177920-02 

-0,6167130-06 

0,4918350-05 

-0,1604270-05 

Z13 

0.1000000*03 

-0.1767900-02 

-0.1025020-02 

-0.3154510-06 

0.0 

-0,6335080-06 
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3.4. 7. 9 Cell Node Strains 


JOB TITU-: U9U TEST CASE - CUBE IN SIMPLE TENSION WITH PLASTICITY 

INTERIOR STRAIN AT TIME = 1.00 


NODE 

EPS >0< 

EPS YY 

EPS Zt 

F''«XY 

EPSXZ 

EPSYZ 

3 

O.lOOOOKHOl 

-0.2500090+00 

-0.2500090+00 

0.0 

-0. 4623090-07 

0.2326870-09 

2 

0.1000010+01 

-0.2500060+00 

-0.2500150+00 

0.0 

-0.346731D-07 

0.174S15D-09 

1 

O.IOCOOIO+QI 

-0.2500090+00 

-0.2500090+00 

0.0 

0.0 

0.0 

11 

0.1000010+01 

-0.2500150+00 

-0.2500050+00 

0.6841120-08 

0.0 

-0.5699870-08 

21 

0.1000010+01 

-0.2500090+00 

-0.2500090+00 

0.1178820-07 

0.0 

-0.7599830-08 

22 

0.1000010+01 

-0.2500060+00 

-0.2500150+00 

0.0 

-0. 2772650-07 

-0.1122520-07 

23 

0.1000010+01 

-0.2500090+00 

-0.2500090+00 

-0.101603D-07 

-0.3696860-07 

-0.736715D-08 

13 

0.1000010+01 

-0.2500150+00 

-0.2500060+00 

-0.7620220-08 

0.0 

-0.5350650-08 

103 

0.1000000+01 

-0.2500060+00 

-0.2500060+00 

0.1838540-08 

-0.2240230-08 

0.0 

lot 

0.1000000+01 

-0.2500060+00 

-0.2500060+00 

-0.4891400-08 

0. 6710600-07 

0.0 

121 

0. 1000000+01 

-0.2500060+00 

-0.2500060+00 

0.1279080-07 

0.7383600-07 

0.0 

123 

0.1000000+01 

* -0.2500060+00 

-0.2500060+00 

-0.1340190-07 

0.1838310-07 

0.0 

203 

0.1000010+01 

-0.2500090+00 

-0.2500090+00 

0.2451390-08 

0.4324390-07 

0.9494920-08 

202 

0.1000010+01 

-0.2500060+00 

-0.2500150+00 

0.0 

0.9953900-07 

0.7121180-08 

201 

0.1000010+01 

-0.2500090+00 

-0.2500090+00 

-0.6521870-08 

0.8947470-07 

0.0 

211 

0.1000010+01 

-0.2500150+00 

-0.2500060+00 

-0. 9416870-09 

0.0 

-0.2216120*07 

221 

0.1000010+01 

-0.2500090+00 

-0.2500090+00 

0.5266290-08 

0.9844800-07 

-0.2954830-07 

222 

0.1000010+01 

-0.2500060+00 

-0.2500150+00 

0.0 

0.1199460-06 

-0.3720130-07 

223 

0. 1000010+01 

-0.2500090+00 

-0.2500090+00 

-0.7708910-08 

0.6147930-07 

-0.2005340-07 

213 

0.1000010+01 

-0.2500150+00 

-0.2500060+00 

-0. 3943130-08 

0.0 

-0. 7918850-08 
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3.4.8 List of Symbols 


Symbol 


Xj 


ui 

«ij 

E 

a 

9 

T 

f 

S 

»o 

X 

yi 

r 

V 

Tijk 





List of Symbols 
Referenced Within Section 3.4 


Description 
Cartesian coordinates 
Inelastic strain rate 
Lame constants 
Displacement rate 
Kronecker delta 
Young's modulus 

Coefficient of thermal expansion 
Time derivative of temperature 
Time derivative of body forces 
Boundary of body to be analyzed 
Point on boundary S 
Integration point 



Interior of body to be analyzed 

Kelvin solution 

Stresses derived from Gi’j 

Tractions derived from T-jjk and 
surface normal 

1/2 6jj if S is smooth at 

otherwise depends on surface geometry at i 


Page 

3.4-6 

3.4-6 

3.4-6 

3.4-6 

3.4-6 

3.4-6 

3.4-6 

3.4-6 

3.4- 6 

3.4- 7 
3.4-7 
3.4-7 
3.4-7 
3.4-7 
3.4-7 
3.4-7 
3.4-7 
3.4-7 

3.4-7 
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List of Symbols 

Referenced Within Section 3.4 (continued) 


Symbol 

Description 

Page 

ni 

Surface normal 

3.4-7 

^i 

Tractions 

3.4-8 

<^jk 

Stress rate 

3.4-8 

Mijk ) 

Higher order kernels derived 
from Gi-j by differentiation 
and use of Hooke's Law 

3.4-8 

— 

Placed over symbol, denotes use of 
a local axis system 

3.4-8 

X 

Vector of all unknown freedoms 
(displacements and tractions) 

3.4-10 

y 

Vector of all known freedoms 
(displacements and tractions) 

3.4-10 

p 

Mass density 

3,4-12 

Cl 

Dilatational wave speed 

3.4-12 

C2 

Distortional wave speed 

3,4-12 

L(-) 

Laplace transform 

3.4-12 

— - 

Denotes (in Section 3. 4. 3. 3) the Laplace 
transform of a function 

3.4-12 

s 

Transform parameter 

3.4-12 

u 

Frequency 

3.4-14 

®ij‘ ^ij 

Laplace transforms of the dynamic 
kernel functions 

3.4-14 

P 

Damping ratio 

3.4-17 

n 

Coefficient of viscosity 

3.4-17 
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